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This  report  includes  a  review  of  previous  work  on  shape  from  shading  and  pho- 
toclinometry.  Novel  features  of  the  new  scheme  are  introduced  one  at  a  time  to  make 
it  easier  to  see  what  each  contributes.  Included  is  a  discussion  of  implementation  de¬ 
tails  that  are  important  if  exact  algebraic  solutions  of  synthetic  shape-from-shading 
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will  lead  to  better  performance  on  real  data. 
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1.  Background 

The  first  method  developed  for  solving  a  shape-from-shading  problem  was 
restricted  to  surfaces  with  special  reflecting  properties  [Rindfleisch  66].  For  the 
surfaces  that  Rindfleisch  considered,  profiles  of  the  solution  can  be  obtained 
by  integrating  along  predetermined  straight  lines  in  the  image.  The  general 
problem  was  formulated  and  solved  later  [Horn  70,  75],  using  the  method 
of  characteristic  strip  expansion  applied  to  the  nonlinear  first-order  partial 
differential  image  irradiance  equation.  When  the  light  sources  and  the  viewer 
are  far  away  from  the  scene  being  viewed,  use  of  the  reflectance  map  makes 
the  analysis  of  shape-from-shading  algorithms  much  easier  [Horn  77]  [Horn 
&;  Sjoberg  79].  Several  iterative  schemes,  mostly  based  on  minimization  of 
some  functional  containing  an  integral  of  the  brightness  error,  arose  later 
[Woodham  77]  [Strat  79]  [Ikeuchi  &  Horn  81]  [Kirk  84,  87]  [Brooks  &  Horn  85] 
[Horn  &  Brooks  86]  [Frankot  &  Chellappa  88]. 

The  new  method  presented  here  was  developed  in  part  as  a  response 
to  recent  attention  to  integr ability1  [Horn  &:  Brooks  86]  [Frankot  &  Chel¬ 
lappa  88]  and  exploits  the  idea  of  a  coupled  system  of  equations  for  depth  and 
slope  [Harris  86,  87]  [Horn  88].  It  borrows  from  well-known  variational  ap¬ 
proaches  to  the  problem  [Ikeuchi  &  Horn  81]  [Brooks  &  Horn  85]  and  an  exist¬ 
ing  least-squares  method  for  estimating  surface  shape  given  a  needle  map  (see 
[Ikeuchi  84],  chapter  11  in  [Horn  86],  and  [Horn  &:  Brooks  86]).  For  one  choice 
of  parameters,  the  new  method  becomes  similar  to  one  of  the  first  iterative 
methods  ever  developed  for  shape  from  shading  on  a  regular  grid  [Strat  79], 
while  it  degenerates  into  another  well-known  method  [Ikeuchi  &  Horn  81]  for  a 
different  choice  of  parameters.  If  the  brightness  error  term  is  dropped,  then  it 
becomes  a  surface  interpolation  method  [Harris  86,  87].  The  computational  ef¬ 
fort  grows  rapidly  with  image  size,  so  the  new  method  can  benefit  from  proper 
multigrid  implementation  [Brandt  77,  80,  84]  [Brandt  &  Dinar  79],  as  can 
existing  iterative  shape-from-shading  schemes  [Terzopolous  83,  84]  [Kirk  84, 
87].  Alternatively,  one  can  apply  so-called  direct  methods  for  solving  Poisson’s 
equations  [Simchony,  Chellappa  &  Shao  89]. 

Linear  expansion  of  the  reflectance  map  about  the  current  • -st’-rate  of 
the  surface  gradient  leads  to  more  rapid  convergence.  More  impor  ly,  this 
modification  often  allows  the  scheme  to  converge  when  simpler  schemes  di¬ 
verge,  or  get  stuck  in  local  minima  of  the  functional.  Most  existing  iterative 
shape-from-shading  methods  handle  only  relatively  simple  surfaces  and  so 
could  benefit  from  a  retrofit  of  this  idea. 

The  new  scheme  was  tested  on  a  number  of  synthetic  images  of  increas¬ 
ing  complexity,  including  some  generated  from  digital  terrain  models  of  steep, 

1 A  gradient  field  is  integrable  if  it  is  the  gradient  of  some  surface  height  function. 


2.  Review  of  Problem  Formulation 


3 


wrinkled  surfaces,  such  as  a  glacial  cirque  with  numerous  gullies.  Shown  in 
Figure  1(a)  is  a  shaded  view  of  a  digital  terrain  model,  with  lighting  from 
the  Northwest.  This  is  the  input  provided  to  the  algorithm.  The  underlying 
231  x  178  digital  terrain  model  was  constructed  from  a  detailed  contour  map, 
shown  in  Figure  2,  of  Huntington  ravine  on  the  eastern  slopes  of  Mount  Wash¬ 
ington  in  the  White  Mountains  of  New  Hampshire2.  Shown  in  Figure  1(b) 
is  a  shaded  view  of  the  same  digital  terrain  model  with  lighting  from  the 
Northeast.  This  is  not  available  to  the  algorithm,  but  is  shown  here  to  make 
apparent  features  of  the  surface  that  may  not  stand  out  as  well  in  the  other 
shaded  view.  Figure  1(c)  shows  a  shaded  view  of  the  surface  reconstructed 
by  the  algorithm,  with  lighting  from  the  Northwest — it  matches  Figure  1(a) 
exactly.  More  importantly,  the  shaded  of  view  of  the  reconstructed  surface 
with  lighting  from  the  Northeast,  shown  in  Figure  1(d),  matches  Figure  1(b) 
exactly3. 

With  proper  boundary  conditions,  the  new  scheme  recovers  surface  ori¬ 
entation  exactly  when  presented  with  noise-free  synthetic  scenes.  Previous 
iterative  schemes  do  not  find  the  exact  solution,  and  in  fact  wander  away 
from  the  correct  solution  when  it  is  used  as  the  initial  guess.  To  obtain  exact 
algebraic  solutions,  several  details  of  the  implementation  have  to  be  care¬ 
fully  thought  through,  as  discussed  in  section  6.  Simple  surfaces  are  easier  to 
process — with  good  results  even  when  several  of  the  implementation  choices 
are  not  made  in  an  optimal  way.  Similarly,  these  details  perhaps  may  be  of 
lesser  importance  for  real  images,  where  other  error  sources  could  dominate. 

In  the  next  few  sections  we  review  image  formation  and  other  elementary 
ideas  underlying  the  usual  formulation  of  the  shape-from-shading  problem. 
Photoclinometry  is  also  briefly  reviewed  for  the  benefit  of  researchers  in  ma¬ 
chine  vision  who  may  not  be  familiar  with  this  field.  We  then  discuss  both 
the  original  and  the  variational  approach  to  the  shape-from-shading  problem. 
Readers  familiar  with  the  basic  concepts  may  wish  to  skip  over  this  material 
and  go  directly  to  section  5,  where  the  new  scheme  is  derived.  For  additional 
details  see  chapters  10  and  11  in  Robot  Vision  [Horn  86]  and  the  collection  of 
papers,  Shape  from  Shading  [Horn  &;  Brooks  89], 

2.  Review  of  Problem  Formulation 

2.1  Image  Projection  and  Image  Irradiance 

For  many  problems  in  machine  vision  it  is  convenient  to  use  a  camera-centered 
coordinate  system  with  the  origin  at  the  center  of  projection  and  the  Z-axis 

2The  gullies  are  steep  enough  to  be  of  interest  to  ice-climbers. 

3For  additional  examples  of  reconstructions  from  shaded  images,  see  section  7. 
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Figure  2.  Contour  map  from  which  the  digital  terrain  model  used 
to  synthesize  Figures  1(a)  and  (b)  was  interpolated.  The  surface  was 
modelled  as  a  thin  plate  constrained  to  pass  through  the  contours  at  the 
specified  elevations.  The  interpolating  surface  was  found  by  solving  the 
biharmonic  equation,  as  described  at  the  end  of  section  5.4. 

aligned  with  the  optical  axis  (the  perpendicular  from  the  center  of  projection 
to  the  image  plane)4.  We  can  align  the  X-  and  K-axes  with  the  image  plane 
x-  and  y-axes.  Let  the  •principal  distance  (that  is,  the  perpendicular  distance 
from  the  center  of  projection  to  the  image  plane)  be  /,  and  let  the  image  plane 
be  reflected  through  the  center  of  projection  so  that  we  avoid  sign  reversal  of 

4In  photoclinometry  it  is  customary  to  use  an  object-centered  coordinate  system. 
This  is  because  surface  shape  can  be  computed  along  profiles  only  when  strong 
additional  constraint  is  provided,  and  such  constraints  are  best  expressed  in  an 
object-centered  coordinate  system.  Working  in  an  object-centered  coordinate 
system,  however,  makes  the  formulation  of  the  shape-from-shading  problem  con¬ 
siderably  more  complex  (see,  for  example,  [Rindfleisch  66]). 
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the  coordinates.  Then  the  perspective  projection  equations  are 

*  =  /y  and  y  =  /  (!) 

The  problem  is  simplified  if  we  assume  that  the  depth  range  is  small  compared 
with  the  distance  of  the  scene  from  the  viewer  (which  is  often  the  case  when 
we  have  a  narrow  field  of  view,  that  is,  when  we  use  a  telephoto  lens).  Then 
we  have 

x^Lx  and  (2) 

for  some  constant  Zo,  so  that  the  projection  is  approximately  orthographic. 
In  this  case  it  is  convenient  to  rescale  the  image  coordinates  so  that  we  can 
write  x  =  X  and  y  =  Y.  For  work  on  shape  from  shading  it  is  also  convenient 
to  use  z,  height  above  some  reference  plane  perpendicular  to  the  optical  axis, 
rather  than  the  distance  measured  along  the  optical  axis  from  the  center  of 
projection. 

If  we  ignore  vignetting  and  other  imaging  system  defects,  then  image  irra- 
diance  E  at  the  point  (x,  y )  is  related  to  scene  radiance  L  at  the  corresponding 
point  in  the  scene  by  [Horn  86] 

E  =  L  j  cos4  a,  (3) 

where  d  is  the  diameter  of  the  lens  aperture,  /  is  the  principal  distance,  and 
the  off-axis  angle  a  is  given  by 


tana  = 


j  y/x2  +  y2. 


(4) 


Accordingly,  image  irradiance5  is  a  multiple  of  the  scene  radiance,  with  the 
factor  of  proportionality  depending  inversely  on  the  square  of  the  /-number6. 
If  we  have  a  narrow  field  of  view,  the  dependence  on  the  off-axis  angle  a 
can  be  neglected.  Alternatively,  we  can  normalize  the  image  by  dividing 
the  observed  image  irradiance  by  cos4  a  (or  whatever  the  actual  vignetting 
function  happens  to  be). 

We  conclude  from  the  above  that  what  we  measure  in  the  image  is  directly 
proportional  to  scene  radiance,  which  in  turn  depends  on  the  strength  and 
distribution  of  illumination  sources,  the  surface  micro- structure  and  surface 
orientation. 

In  order  to  be  able  to  solve  the  shape  from  shading  problem  from  a  single 
image  we  must  assume  that  the  surface  is  uniform  in  its  reflecting  properties. 
If  we  also  assume  that  the  light  sources  are  far  away,  then  the  irradiance  of 


5Grey-levels  are  quantized  estimates  of  image  irradiance. 

6The  /-number  is  the  ratio  of  the  principal  distance  to  the  diameter  of  the  aperture, 
that  is,  f  jd. 
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different  parts  of  the  scene  will  be  approximately  the  same  and  the  incident 
direction  may  be  taken  as  constant.  Finally,  if  we  assume  that  the  viewer 
is  far  away,  then  the  direction  to  the  viewer  will  be  roughly  the  same  for  all 
points  in  the  scene.  Given  the  above,  we  find  that  scene  radiance  does  not 
depend  on  the  position  in  space  of  a  surface  patch,  but  only  on  its  orientation. 

2.2  Specifying  Surface  Orientation 

Methods  for  recovering  shape  from  shading  depend  on  assumptions  about 
the  continuity  of  surface  height  and  its  partial  derivatives.  First  of  all.  since 
shading  depends  only  on  surface  orientation,  we  must  assume  that  the  surface 
is  continuous  and  that  its  first  partial  derivatives  exist.  Most  formulations 
implicitly  also  require  that  the  first  partial  derivatives  be  continuous,  and 
some  even  require  that  second  partial  derivatives  exist.  The  existence  and 
continuity  of  derivatives  lends  a  certain  “smoothness”  to  the  surface  and  allows 
us  to  construct  local  tangent  planes.  We  can  then  talk  about  the  local  surface 
orientation  in  terms  of  the  orientation  of  these  tangent  planes. 

There  are  several  commonly  used  ways  of  specifying  the  orientation  of  a 
planar  surface  patch,  including: 

•  Unit  surface  normal  n  [Horn  h  Brooks  86]; 

•  Point  on  the  Gaussian  sphere  [Horn  84]; 

•  Surface  gradient  ( p,q )  [Horn  77}; 

•  Stereographic  coordinates  (/,  g)  [Ikeuchi  &  Horn  81]; 

•  Dip  and  strike  (as  defined  in  geology)7; 

•  Luminance  longitude  and  latitude8; 

•  Incident  and  emittance  angles  (i  and  e)9; 

For  our  purposes  here,  the  components  of  the  surface  gradient 

dz  j  dz 

p=di  “d  ?=v  (5) 

7Dip  is  the  angle  between  a  given  surface  and  the  horizontal  plane,  while  strike  is 
the  direction  of  the  intersection  of  the  surface  and  the  horizontal  plane.  The  line 
of  intersection  is  perpendicular  to  the  direction  of  steepest  descent. 

8 Luminance  longitude  and  latitude  are  the  longitude  and  latitude  of  a  point  on 
a  sphere  with  the  given  orientation,  measured  in  a  spherical  coordinate  system 
with  the  poles  at  right  angles  to  both  the  direction  towards  the  source  and  the 
direction  toward  the  viewer. 

incident  and  emittance  angles  are  meaningful  quantities  only  when  there  is  a  single 
source;  and  even  then  there  is  a  two-way  ambiguity  in  surface  orientation  unless 
additional  information  is  provided.  The  same  applies  to  luminance  longitude  and 
latitude. 


2.  Review  of  Problem  Formulation 


7 


will  be  most  directly  us<  '  for  specifying  surface  orientation. 

We  can  convert  bt  een  different  representations  easily.  For  example, 
suppose  that  we  are  to  determine  the  unit  surface  normal  given  the  gradient 
components.  We  know  that  if  we  move  a  small  distance  Sx  in  x,  then  the 
change  in  height  is  Sz  =  pSx  (since  p  is  the  slope  of  the  surface  in  the  x 
direction).  Thus  (1,0,  p)T  is  a  tangent  to  the  surface.  If  we  move  a  small 
distance  6y  in  y,  then  the  change  in  height  is  6z  =  q  6y  (since  q  is  the  slope 
of  the  surface  in  the  y  direction).  Thus  (0,  l,q)T  is  also  a  tangent  to  the 
surface.  The  normal  is  perpendicular  to  all  tangents,  thus  parallel  to  the 
cross-product  of  these  particular  tangents,  that  is  parallel  to  (— p,  — y,  1)T. 
Hence  a  unit  normal  can  be  written  in  the  form 

A  =  7; rrb'x  2  (~Pi  ~q' 1)T-  (6) 

y/l  +p2  +  q 2 

Note  that  this  assumes  that  the  2-component  of  the  surface  normal  is  positive. 
This  is  not  really  a  problem  since  we  can  only  see  surface  elements  whose 
normal  vectors  point  within  7r/2  of  the  direction  toward  the  viewer — other 
surface  elements  are  turned  away  from  the  viewer. 

We  can  use  the  same  notation  to  specify  the  direction  to  a  collimated 
light  source  or  a  small  portion  of  an  extended  source.  We  simply  give  the 
orientation  of  a  surface  element  that  lies  perpendicular  to  the  incident  light 
rays.  So  we  can  write10 


s  = 


1 


V1  +P2s+tf 


(-p»,-<hAY 


(7) 


for  some  ps  and  qs. 


2.3  Reflectance  Map 

We  can  show  the  dependence  of  scene  radiance  on  surface  orientation  in  the 
form  of  the  reflectance  map  R(p,q).  The  reflectance  map  can  be  depicted 
graphically  in  gradient  space11  as  a  series  of  nested  contours  of  constant 
brightness  [Horn  77,  86]. 

The  reflectance  map  may  be  determined  experimentally  by  mounting  a 
sample  of  the  surface  on  a  goniometer  stage  and  measuring  its  brightness  un¬ 
der  the  given  illuminating  conditions  for  various  orientations.  Alternatively, 

10There  is  a  small  problem,  however,  with  this  method  for  specifying  the  direction 
toward  the  light  source:  A  source  may  be  “behind”  the  scene,  with  the  direction 
to  the  source  more  than  7r/2  away  from  the  direction  toward  the  viewer.  In  this 
case  the  z-component  of  the  vector  pointing  toward  the  light  source  is  negative. 
11  The  coordinates  of  gradient  space  are  p  and  q,  the  slopes  of  the  surface  in  the  x 
and  y  direction  respectively. 
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one  may  use  the  image  of  a  calibration  object  (such  as  a  sphere)  for  which 
surface  orientation  is  easily  calculated  at  every  point.  Finally,  a  reflectance 
map  may  be  derived  from  a  phenomenological  model,  such  as  that  of  a  Lam¬ 
bertian  surface.  In  this  case  one  can  integrate  the  product  of  the  bidirectional 
reflectance  distribution  function  (BRDF)  and  the  given  distribution  of  source 
brightness  as  a  function  of  incident  angle  [Horn  &  Sjoberg  79]. 

An  ideal  Lambertian  surface  illuminated  by  a  single  point  source  provides 
a  convenient  example  of  a  reflectance  map12.  Here  the  scene  radiance  is  given 
by  R(p,q )  =  ( Eo/n)cosi ,  where  i  is  the  incident  single  (the  angle  between  the 
surface  normal  and  the  direction  toward  the  source),  while  Eq  is  the  irradiance 
from  the  source  on  a  surface  oriented  perpendicular  to  the  incident  rays.  (The 
above  formula  only  applies  when  i  <  7r/2;  the  scene  radiance  is,  of  course,  zero 
for  i  >  ?r/2.)  Now  cosi  =  n  -  s,  so 


R(p,q)  =  ~ 


l  +  PsP  +  qaq _ 

T  -fp2  +$2\/l T  p2, +q* 


as  long  as  the  numerator  is  positive,  otherwise  R(p,  q)  =  0. 


2.4  Image  Irradiance  Equation 

We  are  new  ready  to  write  down  the  image  irradiance  equation 

E(x,y)  =  0  R(p{x,y),q(x,y)),  (9) 

where  E(x,y )  is  the  irradiance  at  the  poi  it  (i ,  y)  in  the  image,  while  R(p,q) 
is  the  radiance  at  the  corresponding  point  in  the  scene,  at  which  p  =  p(x,y) 
and  q  =  q(x ,  y).  The  proportionality  factor  ft  depends  on  the  /-number  of  the 
imaging  system  (and  may  include  a  scaling  factor  that  depends  on  the  units 
in  which  the  instrument  measures  brightness).  It  is  customary  to  rescade 
image  irradiance  so  that  this  proportionality  factor  may  be  dropped.  If  the 
reflectance  map  has  a  unique  global  extremum,  for  example,  then  the  image 
can  be  normalized  in  this  fashion,  provided  that  a  point  can  be  located  that 
has  the  corresponding  surface  orientation13. 

Scene  radiance  also  depends  on  the  irradiance  of  the  scene  and  a  re¬ 
flectance  factor  (loosely  called  albedo  here).  These  factors  of  proportionality 

12Note  that  shape-from-shading  methods  are  most  definitely  not  restricted  to  Lam¬ 
bertian  surfaces.  Such  special  surfaces  merely  provide  a  convenient  pedagogical 
device  for  illustrating  basic  concepts. 

13If  there  is  a  unique  maximum  in  reflected  brightness,  it  is  convenient  to  rescale 
the  measurements  so  that  this  extremum  corresponds  to  E  =  1.  The  same  ap¬ 
plies  when  there  is  a  unique  minimum,  as  is  the  case  for  the  scanning  electron 
microscope  (SEM). 
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can  be  combined  into  one  that  can  be  taken  care  of  by  normalization  of  im¬ 
age  brightness.  Then  only  the  geometric  dependence  of  image  brightness  on 
surface  orientation  remains  in  R(p,  q),  and  we  can  write  the  image  irradiance 
equation  in  the  simple  form 

E(x,y)  =  R(p(x,y),q(x,y))  (10) 

or 

E{x,y)  =  R(zz{x,y),zv{x,y)),  (11) 

where  p  =  zx  and  q  —  zy  are  the  first  partial  derivatives  of  z  with  respect 
to  x  and  y.  This  is  a  first-order  partial  differential  equation  that  is  typically 
nonlinear,  because  the  reflectance  map  in  most  cases  depends  nonlinearly  on 
the  gradient. 


2.5  Reflectance  Map  Linear  in  Gradient 


Viewed  from  a  sufficiently  great  distance,  the  material  in  the  maria  of  the 
moon  has  the  interesting  property  that  its  brightness  depends  only  on  lu¬ 
minance  longitude,  being  independent  of  luminance  latitude  [Hapke  63,  65]. 
When  luminance  longitude  and  latitude  are  related  to  the  incident  and  emit- 
tance  angles,  it  is  found  that  longitude  is  a  function  of  (cos  ij cos  e).  From  the 
above  we  see  that  cost  —  ns,  while  cose  —  n  •  v,  where  v  =  (0,0, 1)T  is  a 
unit  vector  in  the  direction  toward  the  viewer.  Consequently, 

cosi  n-s  1  ,,  , 

-0-+PsP  +  q,q)-  (12) 


cose 


n  •  v 


x/1  +  p2  +  q 2 

Thus  (cost/ cose)  depends  linearly  on  the  gradient  components  p  and  q,  and 
we  can  write 


R(p,q)  =  /(cp  +  *?),  (13) 

for  some  function  /  and  coefficients  c  and  s.  Both  Lommel-Seeliger’s  and 
Hapke’s  functions  fit  this  mold  [Minnaert  61]  [Hapke  63,  65].  (For  a  few  other 
papers  on  the  reflecting  properties  of  surfaces,  see  [Hapke  81,  84]  [Hapke  &: 
Wells  81]  and  the  bibliography  in  [Horn  &  Brooks  89].)  We  can,  without  loss 
of  generality14,  arrange  for  c2  +  s2  =  1. 

If  the  function  /  is  continuous  and  monotonic15,  we  can  find  an  inverse 

cp  +  sq  =  f~l(E(x,y)).  (14) 


14  We  see  that  c  :  s  =  ps  :  q3,  so  that  the  direction  specified  in  the  image  by  (c,s)  is 
the  direction  “toward  the  source,”  that  is,  the  projection  into  the  image  plane  of 
the  vector  s  toward  the  light  source. 

15If  the  function  /  is  not  monotonic,  there  will  be  more  than  one  solution  for 
certain  brightness  values.  In  this  case  one  may  need  to  introduce  assumptions 
about  continuity  of  the  derivatives  in  order  to  decide  which  solution  to  choose. 
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The  slope  in 

the  image  direction  (c,  s)  is 

(15) 

We  can  integrate16  out  this  slope  along  the  line 

*(£)  =  *0  +  and  y(0  =  yo  +  s£, 

(16) 

to  obtain 

z(0  =  Z0  +  J  f  1  y(»7)))  drf. 

(17) 

An  extension  of  the  above  approach  allows  one  to  take  into  account  perspective 
projection  as  well  as  finite  distance  to  the  light  source  [Rindfleisch  66].  Two 
changes  need  to  be  made;  one  is  that  the  reflectance  map  now  is  no  longer 
independent  of  image  position  (since  the  directions  to  the  viewer  and  the 
source  vary  significantly);  and  the  other  that  the  integral  is  for  the  logarithm 
of  the  radial  distance  from  the  center  of  projection,  as  opposed  to  distance 
measured  parallel  to  the  optical  axis. 

The  above  was  the  first  shape-from-shading  or  photoclinometric  problem 
ever  solved  in  other  than  a  heuristic  fashion.  The  original  formulation  was 
considerably  more  complex  than  described  above,  as  the  result  of  the  use  of  full 
perspective  projection,  the  lack  of  the  notion  of  anything  like  the  reflectance 
map,  and  the  use  of  an  object-centered  coordinate  system  [Rindfleisch  66], 

Note  that  we  obtain  profiles  of  the  surface  by  integrating  along  predeter¬ 
mined  straight  lines  in  the  image.  Each  profile  has  its  own  unknown  constant 
of  integration,  so  there  is  a  great  deal  of  ambiguity  in  the  recovery  of  surface 
shape.  In  fact,  if  z(x,y )  is  a  solution,  so  is 

z(x,y)  =  z(x,y)  -hg(sx-cy)  (18) 

for  an  arbitrary  function  gl  This  is  true  because 

=  zx  +  s  g'(s  x  —  cy)  and  zy  =  zy  —  cg'(sx  —  cy),  (19) 
so 

cp  +  sq  =  cp  +  sq,  (20) 

where  p  =  zx  and  q  =  zy.  It  follows  that  R(p,  q)  =  R(p ,  q).  This  ambiguity  can 
be  removed  if  fin  initial  curve  is  given  from  which  the  profiles  can  be  started. 
Such  an  initial  curve  is  typically  not  available  in  practice.  Ambiguity  is  not 
restricted  to  the  special  case  of  a  reflectance  map  that  is  linear  in  the  gradient: 
Without  additional  constta  shape-from-shading  problems  typically  do  not 
have  a  unique  solution. 

16T  he  integration  is,  of  coarse,  .  .  ied  out  numerically,  since  the  integrand  is  derived 
from  image  measurements  and  not  represented  as  an  analytic  function. 
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2.6  Low  Gradient  Terrain  and  Oblique  Illumination 

If  we  are  looking  at  a  surface  where  the  gradient  (p,  q )  is  small,  we  can  ap¬ 
proximate  the  reflectance  map  using  series  expansion 

R{p,  <l)  ^  #(0>  0)  4-  p  Rp( 0, 0)  +  q  Rq( 0, 0).  (21) 

This  approach  does  not  work  when  the  reflectance  map  is  rotationally  sym¬ 
metric,  since  the  first-order  terms  then  drop  out17.  If  the  illumination  is 
oblique,  however,  we  can  apply  the  method  in  the  previous  section  to  get  a 
first  estimate  of  the  surface.  Letting  c  =  i2p(0, 0),  s  =  Rq( 0,0)  and 

r1  (£(*,  y))  =  E(x,  y)  -  R( 0, 0),  (22) 

we  find  that 

z(0  =  *0  +  7  ==  =  f  (E(x(n),yW)  -  R(0,0)\  dr).  (23) 

JK‘(0,0)  +  K‘(0,0)J<>  V  ' 

(For  a  related  frequency  domain  approach  see  [Pentland  88].) 

One  might  imagine  that  the  above  would  provide  a  good  way  to  get  initial 
conditions  for  an  iterative  shape  from  shading  method.  Unfortunately,  this  is 
not  very  helpful,  because  of  the  remaining  ambiguity  in  the  direction  at  right 
angles  to  that  of  profile  integration.  Iterative  methods  already  rapidly  get 
adequate  variations  in  height  along  “down-sun  profiles,”  but  then  struggle  for 
a  long  time  to  try  to  get  these  profiles  tied  together  in  the  direction  at  right 
angles. 

The  above  also  suggests  that  errors  in  gradients  of  a  computed  solution 
axe  likely  to  be  small  in  the  direction  towards  or  away  “from  the  source” 
and  large  in  the  direction  at  right  angles.  It  should  also  be  clear  that  it  is 
relatively  easy  to  find  solutions  for  slowly  undulating  surfaces  (where  p  and  q 
remain  small)  with  oblique  illumination  (as  in  [Kirk  8”]).  It  is  harder  to  deal 
with  cases  where  the  surface  gradient  varies  widely,  and  with  cases  where  the 
source  is  near  the  viewer. 


3.  Brief  Review  of  Photoclinometry 

Photoclinometry  is  the  recovery  of  surface  slopes  from  images  [Wilhelms  64] 
[Rindfleisch  66]  [Lambiotte  &  Taylor  67]  [Watson  68]  [Lucchitta  &  Gambell  70] 
[Tyler,  Simpson  &  Moore  71]  [Rowan,  McCauley  &  Holm  71]  [Bonner  & 
Small  73]  [Wildey  75]  [Squyres  81]  [Howard,  Blasius  &  Cutt  82].  Many  papers 

17The  reflectance  map  is  rotationally  symmetric,  for  example,  when  the  source  is 
where  the  viewer  is,  or  when  an  extended  source  is  symmetrically  distributed 
about  the  direction  toward  the  viewer. 
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and  abstracts  relating  to  this  subject  appear  in  places  that  may  seem  inacces¬ 
sible  to  someone  working  in  machine  vision  [Davis,  Soderblom,  &  Eliason  82] 
[Passey  &  Shoemaker  82]  [Davis  &  McEwen  84]  [Davis  &  Soderblom  83,  84] 
[Malin  &  Danielson  84]  [Wilson  et  al.  84]  [McEwen  85]  [Wilson  et  al.  85]  (For 
additional  references  see  [Horn  h  Brooks  89]).  Superficially,  photoclinometry 
may  appear  to  be  just  another  name  for  shape  from  shading.  Two  different 
groups  of  researchers  independently  tackled  the  problem  of  recovering  surface 
shape  from  spatial  brightness  variations  in  single  images.  Astrogeologists  and 
workers  in  machine  vision  became  aware  of  each  other’s  interests  only  a  few 
years  ago.  The  underlying  goals  of  the  two  groups  are  related,  but  there  are 
some  differences  in  approach  that  may  be  worthy  of  a  brief  discussion. 


3.1  Photoclinometry  versus  Shape  from  Shading 

•  First,  photoclinometry  has  focused  mostly  on  profile  methods  (photo- 
clinometrists  now  refer  to  existing  shape-from-shading  methods  as  area- 
based  photoclinometry,  as  opposed  to  profile -based).  This  came  about 
in  large  part  because  several  of  the  surfaces  of  interest  to  the  astroge- 
ologist  have  reflecting  properties  that  allow  numerical  integration  along 
predetermined  lines  in  the  image,  as  discussed  above  in  section  2.5  [Rind- 
fleisch  66].  Later,  a  similar  profile  integration  approach  was  applied  to 
other  kinds  of  surfaces,  by  using  strong  assumptions  about  local  surface 
geometry  instead.  The  assumption  that  the  surface  is  locally  cylindrical 
leads  to  such  a  profile  integration  scheme  [Wildey  86],  for  example.  More 
commonly,  however,  it  has  been  assumed  that  the  cross-track  slope  is 
zero,  in  a  suitable  object-centered  coordinate  system  [Squyres  81].  This 
may  be  reasonable  when  one  is  considering  a  cross-section  of  a  linearly  ex¬ 
tended  feature,  like  a  ridge,  a  graben,  or  a  central  section  of  a  rotationally 
symmetric  feature  like  a  crater. 

•  The  introduction  of  constraints  that  are  easiest  to  express  in  an  object- 
centered  coordinate  system  leads  away  from  use  of  a  camera-centered 
coordinate  system  and  to  complex  coordinate  transformations  that  tend 
to  obscure  the  underlying  problem.  A  classic  paper  on  photoclinometry 
[Rindfleisch  66]  is  difficult  to  read  for  this  reason,  and  as  a  result  had  little 
impact  on  the  field.  On  the  other  hand,  it  must  be  acknowledged  that 
this  paper  dealt  properly  with  perspective  projection,  which  is  important 
when  the  field  of  view  is  large.  In  all  but  the  earliest  work  on  shape 
from  shading  [Horn  70,  75],  the  assumption  is  made  that  the  projection 
is  approximately  orthographic.  This  simplifies  the  equations  and  allows 
introduction  of  the  reflectance  map. 
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•  The  inherent  ambiguity  of  the  problem  does  not  stand  out  as  obviously 
when  one  works  with  profiles,  as  it  does  when  one  tries  to  fully  reconstruct 
surfaces.  This  is  perhaps  why  workers  on  shape  from  shading  have  been 
more  concerned  with  ambiguity,  and  why  they  have  emphasized  the  im¬ 
portance  of  singular  points  and  occluding  boundaries  [Bruss  82]  [Deift  & 
Sylvester  81]  [Brooks  83]  [Blake,  Zisserman  &  Knowles  85]  [Saxberg  88]. 

•  The  recovery  of  shape  is  more  complex  than  the  computation  of  a  set  of 
profiles.  Consequently  much  of  the  work  in  shape  from  shading  has  been 
restricted  to  simple  shapes.  At  the  same  time,  there  has  been  extensive 
testing  of  shape  from  shading  algorithms  on  synthetic  data.  This  is  some¬ 
thing  that  is  important  for  work  on  shape  from  shading,  but  makes  little 
sense  for  the  study  of  simple  profile  methods,  except  to  test  for  errors  in 
the  procedures  used  for  inverting  the  photometric  function. 

•  Shape-from-shading  methods  easily  deal  with  arbitrary  collections  of  col¬ 
limated  light  sources  and  extended  sources,  since  these  can  be  accommo¬ 
dated  in  the  reflectance  map  by  integrating  the  BRDF  and  the  source 
distribution.  In  astrogeology  there  is  only  one  source  of  light  (if  we  ig¬ 
nore  mutual  illumination  or  interflection  between  surfaces),  so  methods 
for  dealing  with  multiple  sources  or  extended  sources  were  not  developed. 

•  Calibration  objects  are  used  both  in  photoclinometry  and  shape  from 
shading.  In  photoclinometry  the  data  derived  is  used  to  fit  parameters  to 
phenomenological  models  such  as  those  of  Minnaert,  Lommel  and  Seel- 
iger,  Hapke,  and  Lambert.  In  work  on  shape  from  shading  the  numerical 
data  is  at  times  used  directly  without  further  curve  fitting.  The  pa¬ 
rameterized  models  have  the  advantage  that  they  permit  extrapolation 
of  observations  to  situations  not  encountered  on  the  calibration  object. 
This  is  not  an  issue  if  the  calibration  object  contains  surface  elements 
with  all  possible  orientations,  as  it  will  if  it  is  smooth  and  convex. 

•  Normalization  of  brightness  measurements  is  treated  slightly  differently 
too.  If  the  imaging  device  is  linear,  one  is  looking  for  a  single  overall 
scale  factor.  In  photoclinometry  this  factor  is  often  estimated  by  looking 
for  a  region  that  is  more  or  less  flat  and  has  known  orientation  in  the 
object-centered  coordinate  system.  In  shape  from  shading  the  brightness 
of  singular  points  is  often  used  to  normalize  brightness  measurements 
instead.  The  choice  depends  in  part  on  what  is  known  about  the  scene, 
what  the  shapes  of  the  objects  are  (that  is,  are  singular  points  or  oc¬ 
cluding  boundaries  imaged)  and  how  the  surface  reflects  light  (that  is,  is 
there  a  unique  global  extremum  in  brightness). 

•  Finally,  simple  profiling  methods  usually  only  require  continuity  of  the 
surface  and  existence  of  the  first  derivative  (unless  there  is  an  ambiguity 
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in  the  inversion  of  the  photometric  function  whose  resolution  requires  that 
neighboring  derivatives  are  similar).  Most  shape  from  shading  methods 
require  continuous  first  derivatives  and  the  existence  of  second  derivatives 
(In  some  cases  use  is  made  of  the  equality  of  the  second  cross-derivatives 
taken  in  different  order,  that  is,  zxy  =  zyz).  This  means  that  these  meth¬ 
ods  do  not  work  well  on  scenes  composed  of  objecto  that  are  only  piecewise 
smooth,  unless  appropriately  modified18  (but  see  [Malik  &  Maydan  89]). 


3.2  Profiling  Methods 

We  have  seen  in  section  2.5  how  special  photometric  properties  sometimes 
allow  one  to  calculate  a  profile  by  integration  along  predetermined  straight 
lines  in  the  image.  The  other  approach  commonly  used  in  photoclinometry 
to  permit  simple  integration  is  to  make  strong  assumptions  about  the  surface 
shape,  most  commonly  that,  in  a  suitable  object-centered  coordinate  system, 
the  slope  of  the  surface  is  zero  in  a  direction  at  right  angles  to  the  direction  in 
which  the  profile  is  being  computed.  Local  surface  orientation  has  two  degrees 
of  freedom.  The  measured  brightness  provides  one  constraint.  A  second 
constraint  is  needed  to  obtain  a  solution  for  surface  orientation.  A  known 
tangent  of  the  surface  can  provide  the  needed  information.  Two  common 
cases  are  treated  in  astrogeology: 

(a)  features  that  appear  to  be  linearly  extended  (such  as  some  ridges  and 
grabens),  in  a  direction  presumed  to  be  “horizontal”  (that  is,  in  the  local 
average  tangent  plane); 

(b)  features  that  appear  to  be  rotationally  symmetric  (like  craters),  with 
symmetry  axis  presumed  to  be  “vertical”  (that  is,  perpendicular  to  the 
average  local  tangent  plane). 

In  each  case,  the  profile  is  taken  “across”  the  feature,  that  is,  in  a  direction 
perpendicular  to  the  intersection  of  the  surface  with  the  average  local  tangent 
plane.  Equivalently,  it  is  assumed  that  the  cross-track  slope  is  zero  in  the 
object-centered  coordinate  system. 

One  problem  with  this  approach  is  that  we  obtain  a  profile  in  a  plane 
containing  the  viewer  and  the  light  source,  not  a  “vertical”  profile,  one  that  is 
perpendicular  to  the  average  local  tangent  plane.  One  way  to  deal  with  this 
is  to  iteratively  adjust  for  the  image  displacement  resulting  from  fluctuations 
in  height  on  the  surface,  using  first  a  scan  that  really  is  just  a  straight  line  in 

18Methods  for  recovering  the  shapes  of  polyhedral  objects  using  shading  on  the  faces 
and  the  directions  of  the  projections  of  the  edges  into  the  image  are  discussed  in 
[Sugihara  86]  and  [Horn  86]. 
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the  image,  then  using  the  estimated  profile  to  introduce  appropriate  lateral 
displacements  into  the  scan  line,  and  so  on  [Davis  &  Soderblom  84]. 

It  turns  out  that  the  standard  photoclinometric  profile  approach  can 
be  easily  generalized  to  arbitrary  tangent  directions,  ones  that  need  not  be 
perpendicular  to  the  profile,  and  also  to  nonzero  slopes.  All  that  we  need 
to  assume  is  that  the  surface  can  locally  be  approximated  by  a  (general) 
cylinder,  that  is,  a  surface  generated  by  sweeping  a  line,  the  generator,  along 
a  curve  in  space.  Suppose  the  direction  of  the  generator  is  given  by  the  vector 
t  =  (a,  b,  c)T .  Note  that  at  each  point  on  the  surface,  a  line  parallel  to  the 
generator  is  tangent  to  the  surface.  Then,  since  the  normal  is  perpendicular 
to  any  tangent,  we  have  t  •  n  =  0  at  every  point  on  the  surface,  or  just 

ap  +  bq  =  c.  (24) 

This,  together  with  the  equation  E  =  R(p,  q),  constitutes  a  pair  of  equations 
in  the  two  unknowns  p  and  q.  There  may,  however,  be  more  than  one  solution 
(or  perhaps  none)  since  one  of  the  equations  is  nonlinear.  Other  means  must 
be  found  to  remove  possible  ambiguity  arising  from  this  circumstance.  Under 
appropriate  oblique  lighting  conditions,  there  will  usually  only  be  one  solution 
for  most  observed  brightness  values. 

From  the  above  we  conclude  that  we  can  recover  surface  orientation  lo¬ 
cally  if  we  assume  that  the  surface  is  cylindrical,  with  known  direction  of  the 
generator.  We  can  integrate  out  the  resulting  gradient  in  any  direction  we 
please,  not  necessarily  across  the  feature.  Also,  the  generator  need  not  lie 
in  the  local  tangent  plane;  we  can  deal  with  other  situations,  as  long  as  we 
know  the  direction  of  the  generator  in  the  camera- centered  coordinate  system. 
Further  generalizations  are  possible,  since  any  means  of  providing  one  more 
constraint  on  p  and  q  will  do. 

In  machine  vision  too,  some  workers  have  used  strong  local  assumptions 
about  the  surface  to  allow  direct  recovery  of  surface  orientation.  For  example, 
if  the  surface  is  assumed  to  be  locally  spherical,  the  first  two  partial  derivatives 
of  brightness  allow  one  to  recover  the  surface  orientation  [Pentland  84]  [Lee 
&  Rosenfeld  85].  Alternatively,  one  may  assume  that  the  surface  is  locally 
cylindrical  [Wildey  84,  86]  to  resolve  the  ambiguity  present  locally  in  the 
general  case. 

4.  Review  of  Shape  from  Shading  Schemes 

4.1  Characteristic  Strips 

The  original  solution  of  the  general  shape  from  shading  problem  uses  the 
method  of  characteristic  strip  expansion  [Horn  70,  75].  The  basic  idea  is  quite 
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easy  to  explain  using  the  reflectance  map  [Horn  77,  86].  Suppose  that  we  are 
at  a  point  (x,y,  z)T  on  the  surface  and  we  wish  to  extend  the  solution  a  small 
distance  in  some  direction  by  taking  a  step  8x  in  x  and  8y  in  y.  We  need  to 
compute  the  change  in  height  6z.  This  we  can  do  if  we  know  the  components 
of  the  gradient,  because 

8z=p8x  +  q8y.  (25) 

So,  as  we  explore  the  surface,  we  need  to  keep  track  of  p  and  q  in  addition  to 
x,  y  and  z.  This  means  that  we  also  need  to  be  able  to  compute  the  changes 
in  p  and  q  when  we  take  the  step.  This  can  be  done  using 

6p  =  r6x  +  s6y  and  8q  =  s8x  +  t8y ,  (26) 

where  r  =  zxx,  s  =  zxy  —  zyx  and  t  =  zvy  are  the  second  partied  derivatives  of 
the  height.  It  seems  that  we  need  to  now  keep  track  of  the  second  derivatives 
also,  and  in  order  to  do  that  we  need  the  third  partial  derivatives,  and  so  on. 

To  avoid  this  infinite  recurrence,  we  take  another  tack.  Note  that  we 
have  not  yet  used  the  image  irradiance  equation  E(x,y )  =  R(p,q).  To  find 
the  brightness  gradient  we  differentiate  this  equation  with  respect  to  x  and  y 
and  so  obtain 


Ex  —  r  Rp  +  sRq 

and 

Ey  —  3  Rp  + 1  Rq  . 

(27) 

At  this  point  we  exploit  the  fact  that  we 

are  free  to  choose  the  direction  of 

the  step  (8x,6y).  Suppose  that  we 

pick 

8x  =  Rp  8$ 

and 

8y  =  Rq  88, 

(28) 

then,  from  equations  (26)  &:  (27)  we  have 

8p  =  Ex  88, 

and 

8q  =  Ey8(. 

(29) 

This  is  the  whole  “trick.”  We  can  summarize  the  above  in  the  set  of  ordinary 

differential  equations 

x  =  Rp,  y  = 

Rq, 

z=pRp  +  qRq 

*o. 

II 

4  = 

:  Ey, 

(30) 

where  the  dot  denotes  differentiation  with  respect  to  £,  a  parameter  that  varies 
along  a  particular  solution  curve  (the  equations  can  be  rescaled  to  make  this 
parameter  be  arc  length).  Note  that  we  actually  have  more  than  a  mere 
characteristic  curve,  since  we  also  know  the  orientation  of  the  surface  at  all 
points  in  this  curve.  This  is  why  a  particular  solution  is  called  a  characteristic 
strip.  The  projection  of  a  characteristic  curve  into  the  image  plane  is  called 
a  base  characteristic. 

The  base  characteristics  are  predetermined  straight  lines  in  the  image 
only  when  the  ratio  x  :  y  —  Rp  :  Rq  is  fixed,  that  is  when  the  reflectance 
map  is  linear  in  p  and  q.  In  general,  one  cannot  integrate  along  arbitrary 
curves  in  the  image.  Also,  an  initial  curve  is  needed  from  which  to  sprout  the 
characteristics  strips. 
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It  turns  out  that  direct  numerical  implementation  of  the  above  equations 
does  not  yield  particularly  good  results,  since  the  paths  of  the  characteristics 
are  affected  by  noise  in  the  image  brightness  measurements  and  errors  tend 
to  accumulate  along  their  length.  In  particularly  bad  cases,  the  base  charac¬ 
teristics  may  even  cross,  which  does  not  make  any  sense  in  terms  of  surface 
shape.  It  is  possible,  however,  to  grow  characteristic  strips  in  parallel  and  use 
a  so-called  sharpening  process  to  keep  neighboring  characteristics  consistent 
[Horn  70,  75].  This  greatly  improves  the  accuracy  of  the  solution,  since  the 
computation  of  surface  orientation  is  tied  more  closely  to  image  brightness 
itself  rather  than  to  the  brightness  gradient.  This  also  makes  it  possible  to 
interpolate  new  characteristic  strips  when  existing  ones  spread  too  fair  apart, 
and  to  remove  some  when  they  approach  each  other  too  closely. 

4.2  Rotationally  Symmetric  Reflectance  Maps 

One  can  get  some  idea  of  how  the  characteristics  explore  a  surface  by  con¬ 
sidering  the  special  case  of  a  rotationally  reflectance  map,  as  might  apply 
when  the  light  source  is  at  the  viewer  (or  when  dealing  with  scanning  electron 
microscope  (SEM)  images).  Suppose  that 

R(p,q)  ~  /(p2  +  92),  (31) 

then 

Rp  =  2  p/'(p2  +  q 2)  and  Rq  =  2  qf'(p2  +  q 2),  (32) 

and  so  the  directions  in  which  the  base  characteristics  grow  are  given  by 

x  =  kp  and  y  =  kq,  (33) 

for  some  k.  That  is,  in  this  case  the  characteristics  are  curves  of  steepest 
ascent  or  descent  on  the  surface.  The  extrema  of  surface  height  are  sources  and 
sinks  of  characteristics.  These  are  the  points  where  the  surface  has  maxima 
in  brightness. 

This  example  illustrates  the  importance  of  so-called  singular  points.  At 
most  image  points,  as  we  have  seen,  the  gradient  is  not  fully  constrained  by 
image  brightness.  Now  suppose  that  R(p,q )  has  a  unique  global  maximum19, 
that  is 

R{p,q)  <  H(p0,9o)  for  all  (p,q)  jk  (p0,?o).  (34) 

A  singular  point  (xo,yo)  in  the  image  is  a  point  where 

E(x  o ,  yo )  =  R(po ,  qo ).  (35) 

At  such  a  point  we  may  conclude  that  (p, g)  =  ( Po,?o )•  Singular  points  in 
general  are  sources  and  sinks  of  characteristic  curves.  Singular  points  provide 

19The  same  argument  applies  when  the  unique  extremum  is  a  minimum,  as  it  is  in 
the  case  of  scanning  electron  microscope  (SEM)  images. 
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strong  constraint  on  possible  solutions  [Horn  70,  75]  [Bruss  82]  [Brooks  83] 
[Saxberg  88]. 

A  limb  is  the  image  of  an  occluding  boundary ,  where  the  surface  is  tan¬ 
gent  to  the  direction  toward  the  viewer.  It  has  been  suggested  that  occluding 
boundaries  provide  strong  constraint  on  possible  solutions  [Ikeuchi  &  Horn  81] 
[Bruss  82].  As  a  consequence  there  has  been  interest  in  representations  for 
surface  orientation  that  behave  well  near  the  occluding  boundary,  unlike  the 
gradient,  which  becomes  infinite  [Ikeuchi  &  Horn  81]  [Horn  ic  Brooks  86]. 
Recently  there  has  been  some  question  as  to  how  much  constraint  occlud¬ 
ing  boundaries  really  provide,  given  that  singular  points  appear  to  already 
strongly  constrain  the  solution  [Brooks  83]  [Saxberg  88]. 

4.3  Existence  and  Uniqueness 

Amazingly,  questions  of  existence  and  uniqueness  of  solutions  of  the  shape- 
from-shading  problem  have  still  not  been  resolved  satisfactorily.  One  problem 
is  that  it  is  hard  to  say  anything  in  general  that  applies  to  all  possible  re¬ 
flectance  maps.  More  can  be  said  when  specific  reflectance  maps  axe  chosen, 
such  as  ones  that  are  linear  in  the  gradient  or  those  that  are  rotationally 
symmetric  [Bruss  82]. 

It  has  recently  been  shown  that  there  exist  impossible  shaded  images , 
that  is,  images  that  do  not  correspond  to  any  surface  illuminated  in  the  spec¬ 
ified  way  [Szeliski  &:  Horn  89].  It  may  turn  out  that  almost  all  images  with 
multiple  singular  points  are  impossible  in  this  sense  [Saxberg  88].  This  is  an 
important  issue,  because  it  may  help  explain  how  our  visual  system  sometimes 
determines  that  the  surface  being  viewed  cannot  possibly  be  uniform  in  its 
reflecting  properties.  One  can  easily  come  up  with  smoothly  shaded  images, 
for  example,  that  do  not  yield  an  impression  of  shape,  instead  appearing  as 
flat  surfaces  with  spatially  varying  reflectance  or  surface  “albedo.”  (See  also 
Figure  10  in  section  7.2). 

4.4  Variational  Formulations 

As  discussed  above  in  section  2.4,  in  the  case  of  a  surface  with  constant  albedo, 
when  both  the  observer  and  the  light  sources  are  far  away,  surface  radiance 
depends  only  on  surface  orientation  and  not  on  position  in  space  and  the 
image  projection  can  be  considered  to  be  orthographic20.  In  this  case  the 

20The  shape-from-shading  problem  can  be  formulated  and  solved  when  the  viewer 
and  the  light  sources  are  not  at  a  great  distance  [Rindfleisch  66]  [Horn  70,  75], 
but  then  scene  radiance  depends  on  position  as  well  as  surface  orientation,  and 
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image  irradiance  equation  becomes  just 

E(x,y)  =  R(p(x,  y ),  q(x,  y)) ,  (36) 

where  E(x,y)  is  the  image  irradiance  at  the  point  (x,y),  while  R(p,q),  the 
reflectance  map ,  is  the  (normalized)  scene  radiance  of  a  surface  patch  with 
orientation  specified  by  the  partial  derivatives 

dz  j  &Z 

p=Tx  and  «=^,  (37) 

of  surface  height  z(x,y )  above  some  reference  plane  perpendicular  to  the  op¬ 
tical  axis. 

The  task  is  to  find  z(x,y)  given  the  image  E(x,y)  and  the  reflectance 
map  R(p,q).  Additional  constraints,  such  as  boundary  conditions  and  sin¬ 
gular  points,  are  needed  to  ensure  that  there  is  a  unique  solution  [Brass  82] 
[Deift  &:  Sylvester  81]  [Blake,  Zisserman  &  Knowles  85]  [Saxberg  88].  If  we 
ignore  integrability21 ,  some  versions  of  the  problem  of  shape  from  shading 
may  be  considered  to  be  ill-posed22,  that  is,  there  is  not  a  unique  solution 
(p(x,y),<7(x,y)}  that  minimizes  the  brightness  error 

JJ (E(x,y)  -  R(p,q))2  dxdy.  (38) 

In  fact  the  error  can  be  made  equal  to  zero  for  an  infinite  number  of  choices  for 
{p{x,  y),  q(x,  y)}.  We  can  choose  one  of  these  solutions  by  finding  the  one  that 
minimizes  some  functional  such  as  a  measure  of  “departure  from  smoothness” 

JJ(pl+P2y  +  ql+<l2y)dxdy,  (39) 

while  satisfying  the  constraint  E(x,y)  =  R(p,q).  Introducing  a  Lagrange 
multiplier  A(x,y)  to  enforce  the  constraint,  we  find  that  we  have  to  minimize 

JJ {{pi  +  P2y  +  9*  +  gj)  +  A(x,  y)  (E  -  R))  dx  dy.  (40) 

The  Euler  equations  are 

Ap+  \(x,y)Rp  =  0  and  Aq  4-  X(x,y)Rq  =  0.  (41) 

After  elimination  of  the  Lagrange  multiplier  A(x,y),  we  are  left  with  the  pair 
of  equations 

RqAp  =  RpAq  and  E(x,y)  =  R(p,q).  (42) 

There  may  not  be  any  simple  convergent  iterative  scheme  for  this  constrained 
variational  problem  [Horn  &  Brooks  86]  (compare  [Wildey  75]). 

the  notion  of  a  reflectance  map  is  not  directly  applicable. 

21 A  gradient-field  (or  needle-diagram)  {p(x,  y),  g(x,y)}  is  integrable  if  there  exists 
some  surface  height  function  z(x,y)  such  that  p(x,  y)  =  zx(x,y)  and  q(x,y)  = 
Zy(x,  y),  where  the  subscripts  denote  partial  derivatives. 

22If  we  impose  integrability,  and  provide  suitable  boundary  conditions,  then  the 
shape-from-shading  problem  is  definitely  not  ill-posed. 
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We  can  approach  this  problem  another  way  using  the  “departure  from 
smoothness”  measure  in  a  penalty  term  [Ikeuchi  &  Horn  81],  looking  instead 
for  a  minimum  of 

JJ  (iE(x'  y )  “  R(Pi  «)) 2  +  A  (pi  +  p\  +  9r  +  9*))  dx  dy.  (43) 

It  should  be  pointed  out  that  a  solution  of  this  “regularized”  problem  is  not  a 
solution  of  the  original  problem,  although  it  may  be  close  to  some  solution  of 
the  original  problem  [Brooks  85].  In  any  case,  this  variational  problem  leads 
to  the  following  coupled  pair  of  second-order  partial  differential  equations: 

A  Ap  =  ~(E(x,  y )  -  R(p,  q))  Rp(p,  q) 

\Aq  =  -( E(x,y )  -  R(p,q))  Rq(p,q ) 

Using  a  discrete  approximation  of  the  Laplacian  operator23 

{A/}*,  »  £ (7„  -  /«), 

where  /  is  a  local  average24  of  /,  and  e  is  the  spacing  between  picture  cells, 
we  arrive  at  the  set  of  equations 


(44) 

(45) 


K\'pki  =  k\ 'pkl  +  ( E(x,y )  -  R(p,q))Rp(p,q) 
k\'  qkl  =  kX'  qkl  -I-  (E(x ,  y)  -  R(p ,  q))  Rq(p ,  q) 
where  A'  =  A/e2.  This  set  of  equations  suggests  the  iterative  scheme 

p1?+i)  =ri?’  +  ~  *(p<n).9w)).R,(P(’‘,,?,")) 

9i?+1)  =5i?’  +  -  *(P(",.9(")))*.(P(n)>< l*”*) 


(46) 


(47) 


where  the  superscript  denotes  the  iteration  number25. 

From  the  above  it  may  appear  that  R(p,  q),  Rp(p ,  q ),  and  Rq(p ,  5)  should 
be  evaluated  using  the  “old”  values  of  p  and  q.  It  turns  out  that  the  numerical 
stability  of  the  scheme  is  somewhat  enhanced  if  they  are  evaluated  instead  at 
the  local  average  values  p  and  q  [Ikeuchi  &  Horn  81]. 

One  might  hope  that  the  correct  solution  of  the  original  shape-from- 
shading  problem  provides  a  stable  point  for  the  iterative  scheme.  This  is  not 
too  likely,  however,  since  we  are  solving  a  modified  problem  that  includes  a 


23There  are  several  methods  for  approximating  the  Laplacian  operator,  including 
five-point  and  nine-point  approximations.  It  is  well  known  that,  while  the  nine- 
point  approximation  involves  more  computation,  its  lowest-order  error  term  has 
a  higher  order  than  that  of  the  five-point  approximation. 

24Here  k  =  4  when  the  local  average  fki  is  computed  using  the  four  edge-adjacent 
neighbors,  while  k  =  10/3,  when  1/5  of  the  average  of  the  corner-adjacent  neigh¬ 
bors  is  added  to  4/5  of  the  average  of  the  edge-adjacent  neighbors  (see  later). 

25These  equations  are  solved  iteratively  because  the  system  of  equations  is  so  large 
and  because  of  the  fact  that  the  reflectance  map  R{p,q)  is  typically  nonlinear. 
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penalty  term.  Consequently,  an  interesting  question  one  might  ask  about  an 
algorithm  such  as  this,  is  whether  it  will  “walk  away”  from  the  correct  solution 
of  the  original  image  irradiance  equation  E(x,y)  =  R{p,q)  when  this  solution 
is  provided  as  an  initial  condition  [Brooks  85].  The  algorithm  described  here 
does  just  that,  since  it  can  trade  off  a  small  amount  of  brightness  error  against 
an  increase  in  surface  smoothness.  At  the  solution,  we  have  E(x ,  y )  =  J?(p,  q ), 
so  that  the  right  hand  sides  of  the  two  coupled  partial  differential  equations 
(equations  (44))  are  zero.  This  implies  that  if  the  solution  of  the  modified 
problem  is  to  be  equal  the  solution  of  the  original  problem  then  the  Laplacians 
of  p  and  q  must  be  equal  to  zero.  This  is  the  case  for  very  few  surfaces,  just 
those  for  which 

Az(x,y)  =  fc,  (48) 

for  some  constant  k.  While  this  includes  all  harmonic  functions,  it  excludes 
most  real  surfaces,  for  which  adjustments  away  from  the  correct  shape  are 
needed  to  assure  equality  of  the  left  and  right  hand  sides  of  equations  (44) 
describing  the  solution  of  the  modified  problem.  In  general,  this  approach  pro¬ 
duces  solutions  that  are  too  smooth,  with  the  amount  of  distortion  depending 
on  the  choice  of  the  parameter  A.  For  related  reasons,  this  algorithm  does 
well  only  on  simple  smooth  shapes,  and  does  not  perform  well  on  complex, 
wrinkled  surfaces. 

4.5  Recovering  Height  from  Gradient 

In  any  case,  we  are  also  still  faced  with  the  problem  of  dealing  with  the  lack 
of  integrabiliiy,  that  is  the  lack  of  a  surface  z(x,  y)  such  that  p(x,  y )  =  zz(x,  y) 
and  q(x,y)  =  z9(x,y)26.  We  can  try  to  find  the  surface  that  has  partial 
derivatives  zx  and  zy  that  come  closest  to  matching  the  computed  p(x,  y)  and 
q(x,y)  by  minimizing 

JJ ((**  ~  P)2  +(zy  ~  q)2)dxdy.  (49) 

This  leads  to  the  Poisson  equation 

Az  =  pz  4-  qy.  (50) 

Using  the  discrete  approximation  of  the  Laplacian  given  above  (equation  (45)) 
yields 

^ \zki  =  ^ zu  -  ( Pr  +  qy),  (51) 

a  set  of  equations  that  suggests  the  following  iterative  scheme: 

4;+1,=4?’-7((p.  }(»?,+{«,}<.?>)-  <52) 

26The  resulting  gradient  field  is  likely  not  to  be  integrable  because  we  have  not 
enforced  the  condition  py  =  qx,  which  corresponds  to  zxy  =  zyz. 


22 


Height  and  Gradient  from  Shading 


where  the  terms  in  braces  are  numerical  estimates  of  the  indicated  derivatives 
at  the  picture  cell  (£,/). 

The  so-called  natural  boundary  conditions 27  here  axe  just 

czx  +  szy  =cp  +  sq,  (53) 

where  (c,  s)  is  a  normal  to  the  boundary. 

Another  way  of  dealing  with  the  integrability  issue  is  to  directly  minimize 

JJ {(E(x,y)  -  R(p,q))2  +  \(py  -  3x)2)  dxdy.  (54) 

This  leads  to  the  coupled  partial  differential  equations  [Horn  &  Brooks  86] 

A  (pyy  ~  <l*y)  =  (£(*>  v)  ~  R(Pi  ?))  Rp > 

><(qzz  ~  Pyx)  =  {E(x,y)  -  R(p,q))Rq. 

This  set  of  equations  can  also  be  discretized  by  introducing  appropriate  finite 
difference  approximations  for  the  second  partial  derivatives  pyy,  qxx  and  the 
cross  derivatives  of  p  and  q.  An  iterative  scheme  is  suggested  once  one  isolates 
the  center  terms  of  the  discrete  approximations  of  pyy  and  qxx.  This  is  very 
similar  to  the  method  developed  by  Strat,  although  he  arrived  at  his  scheme 
directly  in  the  discrete  domain  [Strat  79].  His  iterative  scheme  avoids  the 
excessive  smoothing  of  the  one  described  early,  but  appears  to  be  less  stable, 
in  the  sense  that  it  diverges  under  a  wider  set  of  circumstances. 

5.  New  Coupled  Height  and  Gradient  Scheme 

The  new  shape-from-shading  scheme  will  be  presented  through  a  series  of 
increasingly  more  robust  variational  methods.  We  start  with  the  simplest, 
which  grows  naturally  out  of  what  was  discussed  in  the  previous  section. 


5.1  Fusing  Height  and  Gradient  Recovery 

One  way  of  fusing  the  recovery  of  gradient  from  shading  with  the  recovery  of 
height  from  gradient,  is  to  represent  both  gradient  (p,  q)  and  height  z  in  one 
variational  scheme  and  to  minimize  the  functional 

J j ((£(*> V)  -  R(p,q))2  +  p({zx  -p)2  +  ( zy  -  g)2))  dxdy.  (56) 

Note  that,  as  far  as  p(x,y)  and  g(x,y),  are  concerned  this  is  an  ordinary  calcu¬ 
lus  problem  (since  no  partial  derivatives  of  p  and  q  appear  in  the  integrand). 
Differentiating  the  integrand  with  respect  to  p(x,y )  and  q(x,y)  and  setting 

27 Natural  boundary  conditions  arise  in  variational  problems  where  no  boundary 
conditions  are  explicitly  imposed  [Courant  &  Hilbert  62]. 
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the  result  equal  to  zero  leads  to 


p  =  zx  4-  ~^(E  -  R)RP, 
q  =  zy  +  -( E  -  R)Rq. 


(57) 


Now  z(x,y )  does  not  occur  directly  in  ( E(x,y )  —  R(p,q))  so  we  actually  just 
need  to  minimize 

JJ((^x-p)2 +(zy~q)2)dxdy,  (58) 

and  we  know  from  the  previous  section  that  the  Euler  equation  for  this  vari¬ 
ational  problem  is  just 

Az=px  +  qv.  (59) 

We  now  have  one  equation  for  each  of  p,  q  and  z. 

These  three  equations  are  clearly  satisfied  when  p  =  zx,  q  =  zy  and 
E  —  R.  That  is,  if  a  solution  of  the  original  shape-from-shading  problem 
exists,  then  it  satisfies  this  system  of  equations  exactly  (which  is  more  than 
can  be  said  for  some  other  systems  of  equations  obtained  using  a  variational 
approach,  as  pointed  out  in  section  4.4).  It  is  instructive  to  substitute  the 
expressions  obtained  for  p  and  q  in  px  +  qy : 

Pz  +  qy  =  +  zw  +  ~  ^(E  —  K)(Rpppx  +  Rpq{py  4-  qx)  4-  Rqqqy )  (60) 

-  ( R2pPx  +  RPRM  +  lx)  +  R\<ly)  +  ( EXRP  +  EyRq)Bigr ). 


Since  A z  =  (pz  +  qy),  we  note  that  the  three  equations  above  axe  satisfied 
when 

( RpPx  +  RpRqiPy  +  ?i)  4-  R%qy)  ~~  (EXRP  +  EyRq)  (61) 

=  ( E  —  R)(RppPX  -|-  Rpq{py  +  qx)  +  Rqqqy). 

This  is  exactly  the  equation  obtained  at  the  end  of  section  4.2  in  [Horn  & 
Brooks  86],  where  an  attempt  was  made  to  directly  impose  integrability  us¬ 
ing  the  constraint  py  =  qx.  It  was  stated  there  that  no  convergent  iterative 
scheme  had  been  found  for  solving  this  complicated  nonlinear  partial  differ¬ 
ential  equation  directly.  The  method  described  here  provides  an  indirect  way 
of  solving  this  equation. 

Note  that  the  natural  boundary  conditions  for  z  are  once  again 

czx  +  szy  =  cp  +  sq,  (62) 

where  (c,  s)  is  a  normal  to  the  boundary. 

The  coupled  system  of  equations  above  for  p,  q  and  z  immediately  sug- 
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gests  an  iterative  scheme 

pir1’ =  {*.}£’  +  ±(e-R)r„ 

AT"  =  <*,}*«*  +  (63! 

r* 

AT"  =  A?  +  j  ({pr)iT"  +  MiT") . 

where  we  have  used  the  discrete  approximation  of  the  Laplacian  for  z  intro¬ 
duced  in  equation  (45).  This  new  iterative  scheme  works  well  when  the  initial 
values  given  for  p,  q  and  z  are  close  to  the  solution.  It  will  converge  to  the  ex¬ 
act  solution  if  it  exists;  that  is,  if  there  exist  a  discrete  set  of  values  {z*/}  such 
that  {pjt/}  and  {?*/}  are  the  discrete  estimate  of  the  first  partial  derivatives 
of  z  with  respect  to  x  and  y  respectively  and  Eki  =  R(pki,qki)- 

In  this  case  the  functional  we  wish  to  minimize  can  actually  be  reduced 
to  zero.  It  should  be  apparent  that  for  this  to  happen,  the  estimator  used  for 
the  Laplacian  must  match  the  sum  of  the  convolution  of  the  discrete  estimator 
of  the  x  derivative  with  itself  and  the  convolution  of  the  discrete  estimator  of 
the  y  derivative  with  itself.  (This  and  related  matters  are  taken  up  again  in 
section  6.2.) 

The  algorithm  can  easily  be  tested  using  synthetic  height  data  Zki-  One 
merely  estimates  the  partial  derivatives  using  suitable  discrete  difference  for¬ 
mulae  and  then  uses  the  resulting  values  pu  and  qki  to  compute  the  synthetic 
image  Eki  =  R{pki,qki)-  This  construction  guarantees  that  there  will  be  an 
exact  solution.  If  a  real  image  is  used,  there  is  no  guarantee  that  there  is  an 
exact  solution,  and  the  algorithm  can  at  best  find  a  good  discrete  approx¬ 
imation  of  the  solution  of  the  underlying  continuous  problem.  In  this  case 
the  functional  will  in  fact  not  be  reduced  exactly  to  zero.  In  some  cases  the 
residue  may  be  large.  This  may  be  the  result  of  aliasing  introduced  when 
sampling  the  image,  as  discussed  in  section  6.5,  or  because  in  fact  the  image 
given  could  not  have  arisen  from  shading  on  a  homogeneous  surface  with  the 
reflectance  properties  and  lighting  as  encoded  in  the  reflectance  map — that 
is,  it  is  an  impossible  shaded  image. 

The  iterative  algorithm  described  in  this  section,  while  simple,  is  not  very 
stable  unless  one  is  close  to  the  exact  solution,  particularly  when  the  surface 
is  complex  and  the  reflectance  map  is  not  close  to  linear  in  the  gradient.  It 
can  be  improved  greatly  by  linearizing  the  reflectance  map.  It  can  also  be 
stabilized  by  adding  a  penalty  term  for  departure  from  smoothness.  This 
allows  one  to  come  close  to  the  correct  solution,  at  which  point  the  penalty 
term  is  removed  in  order  to  prevent  it  from  distorting  the  solution.  We  first 
treat  the  linearization  of  the  reflectance  map. 
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5.2  Linearization  of  Reflectance  Map 


We  can  develop  a  better  scheme  than  the  one  developed  in  the  previous  sec¬ 
tion,  while  preserving  the  apparent  linearity  of  the  equations,  by  approximat¬ 
ing  the  reflectance  map  R(p,  q)  locally  by  a  linear  function  of  p  and  q.  There 
are  several  options  for  choice  of  reference  gradient  for  the  series  expansion,  so 
let  us  keep  it  general  for  now  at  (po, go)28-  We  have 

R(p,q)  «  R(po ,  qo )  +  (p  ~  Po)  RP(po,q0)  +  (?  ~  qo)Rq(po,qo)  +  •  •  •  (64) 

Again,  gathering  all  of  the  terms  in  p*/  and  qki  on  the  left  hand  sides  of  the 
equations,  we  now  obtain 


(p  +  R2P) pu  +  RpRq  qki  =  pzx  +  (E  -  R  -  poRP  -  qoRg)RP, 
RqRpPkl  +  (p  +  R\)qkl  =  pzy  +(E  -  R-  poRp  -  q0Rq)Rq, 


(65) 


while  the  equation  for  z  remains  unchanged.  (Note  that  here  R,  Rp  and  Rq 
denote  quantities  evaluated  at  the  reference  gradient  (po,  qo))- 

It  is  convenient  to  rewrite  these  equations  in  terms  of  quantities  relative 
to  the  reference  gradient  (po,qo)-  Let 

6pu  =  Pkl  -  Po  and  6qki  =  qkl  ~  qo, 

6zx  =  zx  -  po  and  8zy  =  zy  —  q0. 

This  leads  to 

(p  +  R2p)  6Pki  +  RpRq  6qki  —  p8zx  +  (E  ~  R)RP, 

RpRq  Sqki  +  (p  +  R])8qu  =  P  8zy  -f  (E  -  R)Rq. 

(The  equations  clearly  simplify  somewhat  if  we  choose  (zx,  zy)  for  the  reference 
gradient  (po,  qo)-)  We  can  view  the  above  as  a  pair  of  linear  equations  for  Spki 
and  8qki.  The  determinant  of  the  2x2  coefficient  matrix, 

D  =^  +  Rl  +  R]),  (68) 

is  always  positive,  so  there  is  no  problem  with  singularities.  The  solution  is 
given  by 


(66) 


(67) 


D  Spki  =  (p  +  R*)  A  -  RpRq  B, 
D  Sqki  =  (/j  +  Rl)  B  -  RqRp  A, 


(69) 


where 


(70) 


A  =  p6zx  +  (E  -  R)RP, 

B  =  p  6zy  +  (E  —  R)Rq. 

This  leads  to  a  convenient  iterative  scheme  where  the  new  values  are  given  by 

p[ni+1)  =  Pon)  +  tplT  and  ?u+1)  =  q(on)  +  tiqff’  (71) 


28The  reference  gradient  will,  of  course,  be  different  at  every  picture  cell,  but  to 
avoid  having  subscripts  on  the  subscripts,  we  will  simple  denote  the  reference 
gradient  at  a  particular  picture  cell  by  (po.tfo)- 
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in  terms  of  the  old  reference  gradient  and  the  increments  computed  above. 
The  new  version  of  the  iterative  scheme  does  not  require  a  great  deal  more 
computation  than  the  simpler  scheme  introduced  in  section  4.5,  since  the 
partial  derivatives  Rp  and  Rq  are  needed  in  any  case. 

5.3  Incorporating  Departure  from  Smoothness  Term 

We  now  introduce  a  penalty  term  for  departure  from  smoothness,  effectively 
combining  the  iterative  method  of  [Ikeuchi  &  Horn  81]  for  recovering  p  and 
q  from  E(x,y)  and  R(p,q),  with  the  scheme  for  recovering  z  given  p  and  q 
discussed  in  section  4.5.  (For  the  moment  we  do  not  linearize  the  reflectance 
map;  this  will  be  addressed  in  section  5.6).  We  look  directly  for  a  minimum 
of 

JJ((E{x,y)-R(p,,))1 

+  M  pi  +  pi  +  +  ll)  (72) 

+  p((zx-p)2  +  (zy  -q)2)')  dxdy. 

The  Euler  equations  of  this  calculus  of  variations  problem  lead  to  the  following 
coupled  system  of  second-order  partial  differential  equations: 

A  Ap  =  ~(E  -  R)RP  -  p(zt  -  p), 

A  A q  — -(E  -  R)Rq  -  p(zy  —  q),  (73) 


Az  =  px  -(-  qy. 

A  discrete  approximation  of  these  equations  can  be  obtained  by  using  the 
discrete  approximation  of  the  Laplacian  operator  introduced  in  equation  (45): 

^  ( Pki  ~  Pki )  =  ~{E  -  R)RP  -  p(zr  -  Pk, ), 

(Qk!  ~  qki)  =  ~(E  -  R)Rq  -  n(zy  -  qkl ),  (74) 

~2  ~  Zkl)  =Pr  +9y 

where  E,  R,  Rp,  and  Rq  are  the  corresponding  values  at  the  picture  cell  ( k ,  l), 
while  zx,  zy,  pz  and  qy  are  discrete  estimates  of  the  partial  derivative  of  z,  p 
and  q  there.  We  can  collect  the  terms  in  pki,  qkl  and  zki  on  one  side  to  obtain 
(kA'  -I-  p)pki  =  {nX'pki  +  pzz)  +  (E  -  R)RP, 

(« A'  +  p)  qki  =  (kX  qkl  +  pzy)  +  (E  -  R)Rq,  (75) 

K>  /C 

Zk{  ~  (pr  +  qy), 

where  A'  =  A/e2.  These  equations  immediately  suggest  an  iterative  scheme, 
where  the  right  hand  sides  are  computed  using  the  current  values  of  the  zku 
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Pkii  and  q/si ,  with  the  results  then  used  to  supply  new  values  for  the  unknowns 
appearing  on  the  left  hand  sides 

From  the  above  it  may  appear  that  R(p,q),  Rp(p,q),  and  Rq(p,q)  should 
be  evaluated  using  the  “old”  values  of  p  and  q.  One  might,  on  the  other  hand, 
argue  that  the  local  average  values  p  and  <j,  or  perhaps  even  the  gradient 
estimates  zx  and  zy,  are  more  appropriate.  Experimentation  suggests  that 
the  scheme  is  most  stable  when  the  local  averages  p  and  q  are  used. 

The  above  scheme  contains  a  penalty  term  for  departure  from  smooth¬ 
ness,  so  it  may  appear  that  it  cannot  to  converge  to  the  exact  solution.  Indeed, 
it  appears  as  if  the  iterative  scheme  will  “walk  away”  from  the  correct  solution 
when  it  is  presented  with  it  as  initial  conditions,  as  discussed  in  section  4.4.  It 
turns  out,  however,  that  the  penalty  term  is  needed  only  to  prevent  instability 
when  far  from  the  solution.  When  we  come  closer  to  the  solution,  A'  can  be 
reduced  to  zero,  and  so  the  penalty  term  drops  out.  It  is  tempring  to  leave 
the  penalty  term  out  right  from  the  start,  since  this  simplifies  the  equations 
a  great  deal,  as  shown  in  section  5.1.  The  contribution  from  the  penalty  term 
does,  however,  help  damp  out  instabilities  when  far  from  the  solution  and  so 
should  be  included.  This  is  particularly  important  with  real  data,  where  one 
cannot  expect  to  find  an  exact  solution. 

(Note  also  that  the  coupled  second  order  partial  differential  equations 
above  (equation  (76))  me  eminently  suited  for  solution  by  coupled  resistive 
grids  [Horn  88].) 


5.4  Relationship  to  Existing  Techniques 

•  Recently  a  new  method  has  been  developed  that  combines  the  itera¬ 
tive  scheme  discussed  in  section  4.4  for  recovering  surface  orientation 
from  shading  with  a  projection  onto  the  subspace  of  integrable  gradients 
[Frankot  &  Chellappa  88).  The  approach  there  is  to  alternately  take  one 
step  of  the  iterative  scheme  [Ikeuchi  &:  Horn  81]  and  to  find  the  “nearest” 
integrable  gradient.  This  gradient  is  then  provided  as  initial  conditions 
for  the  next  step  of  the  iterative  scheme,  ensuring  that  the  gradient  field 
never  departs  far  from  integrability.  The  integrable  gradient  closest  to  a 
given  gradient  field  is  found  using  orthonormal  series  expansion  and  by 
exploiting  the  fact  that  differentiation  in  the  spatial  domain  corresponds 
to  multiplication  by  frequency  in  the  transform  domain. 

•  Similar  results  can  be  obtained  by  using  instead  the  method  described 
in  section  4.5  for  recovering  the  height  z(x,  y )  that  best  matches  a  given 
gradient.  The  resulting  surface  can  then  be  (numerically)  differentiated 
to  obtain  initial  values  for  p{x,y)  and  q(x,y)  for  the  next  step  of  the 
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iterative  scheme.  This  works,  but  not  as  well  as  the  new  scheme  described 
in  the  previous  section. 

•  Next,  note  that  we  obtain  the  scheme  of  [Ikeuchi  &  Horn  81]  (who  ignore 
the  integrability  problem)  discussed  in  section  4.4,  if  we  drop  the  depar¬ 
ture  from  integrability  term  in  the  integrand — that  is,  when  p  =  0.  If  we 
instead  remove  the  departure  from  smoothness  term  in  the  integrand — 
that  is,  when  A  =  0 — we  obtain  something  reminiscent  of  the  iterative 
scheme  of  [Strat  79],  although  Strat  dealt  with  the  integrability  issue  in 
a  different  way. 

•  Finally,  if  we  drop  the  brightness  error  term  in  the  integrand,  we  obtain 
the  scheme  of  [Harris  86,  87]  for  interpolating  from  depth  and  slope.  He 
minimizes 

Jj  (A  (Pi  +  P?  +  9*  +  )  +  ((**  ~  Pf  +  ( 2v  -  ?)2))  dx  dV ■  (76) 

and  arrives  at  the  Euler  equations 

A  A p  =  ~(zx  -p), 

A  A q=-(Zy-q),  (77) 

A  z=px+qy. 

Now  consider  that 

A(Az)  =  A(px  -f  qy).  (78) 

Since  application  of  the  Laplacian  operator  and  differentiation  commute 
we  have 

A(Az)  =  (Ap)r  +  (A?)„  (79) 

or 

A  A(Az)  =  —{ztz  —  Px)  (z yy  —  5y),  (80) 

and  so 

A  A(Az)  =  -Az  +  (px  +  qy)  =  0.  (81) 

So  his  method  solves  the  biharmonic  equation  for  z,  by  solving  a  coupled 
set  of  second-order  partial  differential  equations.  It  does  it  in  an  elegant, 
stable  way  that  permits  introduction  of  constraints  on  both  height  z  and 
gradient  (p,  q).  This  is  a  good  method  for  interpolating  from  sparse  depth 
and  surface  orientation  data. 

The  biharmonic  equation  has  been  employed  to  interpolate  digital  terrain 
models  (DTMs)  from  contour  maps.  Such  DTMs  were  used,  for  example, 
in  [Horn  &  Bachman  78]  [Horn  81]  [Sjoberg  &  Horn  83].  The  obvious  im¬ 
plementations  of  finite  difference  approximations  of  the  biharmonic  operator, 
however,  tend  to  be  unstable  because  some  of  the  weights  are  negative,  and 
because  the  corresponding  coefficient  matrix  lacks  diagonal  dominance.  Also, 
the  treatment  of  boundary  conditions  is  complicated  by  the  fact  that  the 
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support  of  the  biharmonic  operator  is  so  large.  The  scheme  described  above 
circumvents  both  of  these  difficulties — it  was  used  to  interpolate  the  digital 
terrain  model  used  for  the  example  illustrated  by  Figure  1. 

5.5  Boundary  Conditions  &  Nonlinearity  of  Reflectance  Map 

So  far  we  have  assumed  that  suitable  boundary  conditions  axe  available,  that 
is,  the  gradient  is  known  on  the  boundary  of  the  image  region  to  which  the 
computation  is  to  be  applied.  If  this  is  not  the  case,  the  solution  is  likely  not 
to  be  unique.  We  may  nevertheless  try  to  find  a  solution  by  imposing  so-called 
natural  boundary  conditions  [Courant  h  Hilbert  62].  The  natural  boundary 
conditions  for  the  variational  problem  described  here  can  be  shown  to  be 

cpx+spy  =  0  and  cqx+sqs  =  0  (82) 

and 

czx  +  szy  =cp  +  sq  (83) 

where  (c,  s)  is  a  normal  to  the  boundary.  That  is,  the  normal  derivative  of 
the  gradient  is  zero  and  the  normal  derivative  of  the  height  has  to  match  the 
slope  in  the  normal  direction  computed  from  the  gradient. 

In  the  above  we  have  approximated  the  original  partial  differential  equa¬ 
tions  by  a  set  of  discrete  equations,  three  for  every  picture  cell  (one  each  for 
p,  q  and  z ).  If  these  equations  were  linear,  we  could  directly  apply  all  the 
existing  theory  relating  to  convergence  of  various  iterative  schemes  and  how 
one  solves  such  equations  efficiently,  given  that  the  corresponding  coefficient 
matrices  axe  sparse29.  Unfortunately,  the  equations  are  in  general  not  linear, 
because  of  the  nonlinear  dependence  of  the  reflectance  map  R(p ,  q)  on  the  gra¬ 
dient.  In  fact,  in  deriving  the  above  simple  iterative  scheme,  we  have  treated 
f2(p,q),  and  its  derivatives,  as  constant  (independent  of  p  and  q)  during  any 
particular  iterative  adjustment  of  p  and  q. 


5.6  Local  Linear  Approximation  of  Reflectance  Map 

In  section  5.2  we  linearized  the  reflectance  map  in  order  to  improve  the  simple 
scheme  developed  in  section  5.1.  We  now  do  the  same  for  the  more  complex 
scheme  described  in  section  5.4.  We  have 

R(p,q)  «  R(po,qo)  +  (p  -  Po)  Rp(po,qo)  +  (q  ~  qo)  Rg(po,qo)  +  •  ■  •  (84) 

Again,  gathering  all  of  the  terms  in  pki  and  qki  on  the  left  hand  sides  of  the 
equations,  we  now  obtain 

(A"  +  R?p)pki  +  RpRq  qki 

29See  [Lee  88]  for  a  proof  of  convergence  of  an  iterative  shape-from-shading  scheme. 
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=  (K^'Pki  +  pzx)  +  (E  -  R-  p0Rp  -  q0Rq)Rp , 

RqRp  Pki  +  (A"  +  R2q)  qu  (85) 

=  (K^Qkt  +  pzy)  +  (E  -  R-  poRp  -  q0Rq)Rg, 
while  the  equation  for  z  remains  unchanged.  (Note  that  here  R,  Rp  and  Rq 
again  denote  quantities  evaluated  at  the  reference  gradient  (po,  go))-  In  the 
above  we  have  abbreviated  A"  =  kX'  -f  p. 

It  is  convenient  to  rewrite  these  equations  in  terms  of  quantities  defined 
relative  to  the  reference  gradient: 


8pki  —  Pkt  -  Po 

and 

8qkt  =  -  qo 

tpki  =  Pkt  ~  Po 

and 

&lkl  —  Qki  ~~  Qo 

8zx  =  zx  -  po 

and 

8zy  Zy  go 

This  yields 


(86) 


(A"  -f  Rp)  Spu  +  RpRq  &qki  =  kX'  6pkl  +  p  Szx  +  {E  —  R)RP, 

RpRq  8qki  +  (A"  +  R2)  6qki  —  kX'  8qki  +  p  8zy  +  (E  —  R)Rq. 

(The  equations  clearly  simplify  somewhat  if  we  choose  either  p  and  q  or  zx 
and  zy  for  the  reference  gradient  p<j  and  go.)  We  can  view  the  above  as  a  pair 
of  linear  equations  for  8pki  and  8qki-  The  determinant  of  the  2x2  coefficient 
matrix 

£>  =  \"(\"  +  R2  +  R2)  (88) 

is  always  positive,  so  there  is  no  problem  with  singularities.  The  solution  is 
given  by 

D  8pkl  =  (X"  +  Rq)  A  -  RpRq  B, 


D  8qki  =  (A"  +  R2)  B  —  RqRp  A, 


(89) 


where 


(90) 


A  —  kX'  8pkl  +  p  8zx  +  ( E  —  R)RP, 

B  =  kX'  8qkl  +  p8zy  +  (E  -  R)Rq. 

This  leads  to  a  convenient  iterative  scheme  where  the  new  values  are  given  by 

Pkl  ~Po  +  SPki  ^  Qki  -Qo  +  SQkl  »  (91) 

in  terms  of  the  old  reference  gradient  and  the  increments  computed  above.  It 
has  been  determined  empirically  that  this  scheme  converges  under  a  far  wider 
set  of  circumstances  than  the  one  presented  in  the  previous  section. 

Experimentation  with  different  reference  gradients,  including  the  old  val¬ 
ues  of  p  and  q ,  the  local  average  p  and  q,  as  well  as  zx  and  zy  showed  that 
the  accuracy  of  the  solution  and  the  convergence  is  affected  by  this  choice.  It 
became  apparent  that  if  we  do  not  want  the  scheme  to  “walk  away”  from  the 
correct  solution,  then  we  should  use  the  old  value  of  p  and  q  for  the  reference 
po  and  go- 
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6.  Some  Implementation  Details 

6.1  Derivative  Estimators  and  Staggered  Grids 

In  one  dimension,  it  is  well-known  from  numerical  analysis  that  the  best  finite 
difference  estimators  of  even  derivatives  have  odd  support,  while  the  best 
estimators  of  odd  derivatives  have  even  support.  These  estimators  are  “best” 
in  the  sense  that  their  lowest-order  error  terms  have  a  small  coefficient  and 
that  they  do  not  attenuate  the  higher  frequencies  as  much  as  the  alternative 
ones.  A  good  estimator  of  the  second  derivative  of  z,  for  example,  is 

{zxx}k  «  \{zk-x  ~  2 zk  +  z*+i),  (92) 

while  a  good  estimator  of  the  first  derivative  of  z  is  just 

{zz}k  «  -(z/fe+i  -  Zk)-  (93) 

Note  that  the  latter,  like  other  estimators  with  even  support  for  odd  deriva¬ 
tives,  gives  an  estimate  valid  at  the  point  midway  between  samples. 

This  suggests  that  one  should  use  staggered  grids.  That  is,  the  arrays 
containing  sampled  values  of  p  and  q  (and  hence  image  brightness  E)  should 
be  offset  by  1/2  pixel  in  both  X  and  y  from  those  for  z  (see  Figure  3).  This 
also  means  that  if  the  image  is  rectangular  and  contains  n  x  m  pixels,  then 
the  array  of  heights  should  be  of  size  (n  +  l)x(m  +  l).  Appropriate  two- 
dimensional  estimators  for  the  first  partial  derivatives  of  z  then  axe  (see  also 
[Horn  &  Schunck  81]): 

{Zr}k,j  K  7r(zk,l+ 1  -  zk,l  +  Zk+\,l+l  -  Zk+l,l) 

1  (94) 

{*»}*, J  W  2l(zk+l,l  -  zk,l  +  zk+l,l+l  -  Zk,l+ 1) 

These  can  be  conveniently  shown  graphically  in  the  form  of  the  stencils: 


The  results  obtained  apply  to  the  point  (fc  +  1/2,  f-f  1/2)  in  the  grid  of  discrete 
values  of  z;  or  the  point  ( k ,  l)  in  the  (offset)  discrete  grid  of  values  of  p  and  q. 
Similar  schemes  can  be  developed  for  the  first  partial  derivatives  of  p  and  q 
needed  in  the  algorithms  introduced  here,  with  the  offsets  now  acting  in  the 
opposite  direction. 
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Figure  3.  It  is  convenient  to  have  the  discrete  grid  for  p,  q  (and  hence 
for  the  image  E  itself)  offset  by  1/2  pixel  in  x  and  1/2  pixel  in  y  from  the 
grid  for  z. 


6.2  Discrete  Estimators  of  the  Laplacian 


We  also  need  to  obtain  local  averages  based  on  discrete  approximations  of  the 
Laplacian  operators.  We  could  simply  use  one  of  the  stencils 


The  second,  diagonal,  form  has  a  higher  coefficient  on  the  lowest-order  error 
term  than  the  first,  edge-adjacent  form,  and  so  is  usually  not  used  by  itself. 
The  diagonal  form  is  also  typically  not  favored  in  iterative  schemes  for  solving 
Poisson’s  equations,  since  it  does  not  suppress  certain  high  frequency  compo¬ 
nents.  We  can  write  a  stencil  for  a  linear  combination  of  the  edge-adjacent 
and  the  diagonal  versions  in  the  form 


a 

4 

1— a 

4 

a 

4 

4 

1  —a 

1 

1— a 

(a  +  1)  e2 

4 

4 

a 

1  —  a 

a 

4 

4 

4 

A  judiciously  chosen  weighted  average,  namely  one  for  which  a  =  1/5,  is 
normally  preferred,  since  this  combination  cancels  the  lowest-order  error  term. 
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If  we  wish  to  prevent  the  iterative  scheme  from  “walking  away”  from  the 
solution,  however,  we  need  to  make  our  estimate  of  the  Laplacian  consistent 
with  repeated  application  of  our  estimators  for  the  first  partial  derivatives. 
That  is,  we  want  our  discrete  estimate  of  A z  to  be  as  close  as  possible  to  our 
discrete  estimate  of 

(zx)x  +  ( zy)y  (95) 

It  is  easy  to  see  that  the  sum  of  the  convolution  of  the  discrete  estimator  for 
the  x-derivative  with  itself  and  the  convolution  of  the  discrete  estimator  for 
the  y-derivative  with  itself  yields  the  diagonal  pattern.  So,  while  the  diagonal 
pattern  is  usually  not  favored  because  it  leads  to  less  stable  iterative  schemes, 
it  appears  to  be  desirable  here  to  avoid  inconsistencies  between  discrete  es¬ 
timators  of  the  first  and  second  partial  derivatives.  Experimentation  with 
various  linear  combinations  bears  this  out.  The  edge-adjacent  stencil  is  very 
stable  and  permits  over-relaxation  (SOR)  with  a  =  2  (see  next  section),  but 
leads  to  some  errors  in  the  solution  with  noisefree  input  data.  The  diagonal 
form  is  less  stable  and  requires  a  reduced  value  for  a,  but  allows  the  scheme  to 
converge  to  the  exact  algebraic  solution  to  problems  that  have  exact  solutions. 

The  incipient  instability  inherent  in  use  of  the  diagonal  form  is  a  reflection 
of  the  fact  that  if  we  think  of  the  discrete  grid  as  a  checkerboard,  then  the 
“red”  and  the  “black”  squares  are  decoupled30.  That  is,  updates  of  red  squares 
are  based  only  on  existing  values  on  red  squares,  while  updates  of  black  squares 
are  based  only  on  existing  values  on  black  squares.  Equivalently,  note  that 
there  is  no  change  in  the  incremental  update  equations  when  we  add  a  discrete 
function  of  the  form 

Szu  =  (-1)*+'  (96) 

to  the  current  values  of  the  height.  The  reason  is  that  the  estimators  of  the 
first  derivatives  and  the  diagonal  form  of  the  Laplacian  estimator  are  com¬ 
pletely  insensitive  to  components  of  this  specific  (high)  spatial  frequency31. 
Fortunately,  the  iterative  update  cannot  inject  components  of  this  frequency 
either,  so  that  if  the  average  of  the  values  of  the  “red”  cells  initially  matches 
the  average  of  the  values  of  the  “black”  cells,  then  it  will  continue  to  do  so. 
The  above  has  not  turned  out  to  be  an  important  issue,  since  the  iteration 

30The  “red”  and  “black”  squares  are  the  cells  for  which  the  sum  of  the  row  and 
column  indexes  are  even  and  odd  respectively. 

31  It  may  appear  that  this  difficulty  stems  from  the  use  of  staggered  grids.  The 
problem  is  even  worse  when  aligned  grids  are  used,  however,  because  the  discrete 
estimator  of  the  Laplacian  consistent  with  simple  central  difference  estimators  of 
the  first  partial  derivatives  has  a  support  that  includes  only  cells  that  are  2e  away 
from  the  center.  And  this  form  of  the  Laplacian  operator  is  known  to  be  badly 
behaved.  We  find  that  there  are  four  decoupled  subsets  of  cells  in  this  case. 
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appears  to  be  stable  with  the  diagonal  form  of  the  average,  that  is,  for  a  =  1, 
when  the  natural  boundary  conditions  are  implemented  with  care. 


6.3  Boundary  Conditions 


The  boundary  conditions  have  also  to  be  dealt  with  properly  to  assure  con¬ 
sistency  between  first-  and  second-order  derivative  estimators.  In  a  simple 
rectangular  image  region,  the  natural  boundary  conditions  for  z  could  be 
implemented  by  simply  taking  the  average  of  the  two  nearest  values  of  the 
appropriate  gradient  component  and  multiplying  by  e  to  obtain  an  offset  from 
the  nearest  value  of  z  in  the  interior  of  the  grid.  That  is,  for  1  <  k  <  n  and 
1  <  1  <  m,  we  could  use 


zk,  0  —  zk,l  —  2^k~ 1,0  Pm) 

zk,m  ~  zk,m  —  1  "b  r  (pit  — 1  ,m  — 1  "b  Pk,m  —  l) 

e  (97) 

Z0,l  =  Z\ ,l  -  +  ?o ,/) 

zn,l  —  Zn-iJ  +  ~(qn-l,l~l  +  <7n-I ,/) 

on  the  left,  right,  bottom  and  top  border  of  a  rectangular  image  region  (the 
corners  are  extrapolated  diagonally  from  the  nearest  point  in  the  interior  using 
both  components  of  the  gradient).  But  this  introduces  a  connection  between 
the  “red” and  the  “black”  cells,  and  so  must  be  in  conflict  with  the  underlying 
discrete  estimators  of  the  derivatives  that  are  being  used. 

One  can  do  better  using  offsets  from  cells  in  the  interior  that  lie  in  diag¬ 
onal  direction  from  the  ones  on  the  boundary.  That  is,  for  2  <  k  <  n  —  1  and 
2  <  /  <  m  —  1,  we  use 


Zk, o  =  -(zfc-1,1  -  e(pk-i,0  ~  qk- 1,0 )  +  Zk+i,i  -  c(pm  +  qk, o)) 

Zk,m  =  2  (2fc-l,m-l  +  e(pjfc_i|Tn_i  +  qk- l,m  — l)  +  2jfc+l,m-l  +  t(pk,m- 1  —  Qk,m-l 

zoj  =  -  (zi.j-i  +  c(po,«-i  —  qo,i-\)  +  zi,/+i  —  «(po,j  +  qo,i)) 

zn,l  =  ^  (zn-l,l-l  +  e(Prt  — 1 1  +  ^n  — 1 ,/  — 1 )  +  Zn-l,l+l  ~  e(pn-l,l  ~  ?n-l,/)) 


(98) 

on  the  left,  right,  bottom  and  top  border  of  a  rectangular  image  region.  The 
corners  are  again  extrapolated  diagonally  from  the  nearest  point  in  the  interior 
using  both  components  of  the  gradient.  Note  that,  in  this  scheme,  one  point 
on  each  side  of  the  corner  has  to  be  similarly  interpolated,  because  only  one 
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of  the  two  values  needed  by  the  above  diagonal  template  lies  in  the  interior 
of  the  region. 

If  the  surface  gradient  is  not  given  on  the  image  boundary,  then  natural 
boundary  conditions  must  be  used  for  p  and  q  as  well.  The  natural  boundary 
condition  is  that  the  normal  derivatives  of  p  and  q  are  zero.  The  simplest 
implementation  is  perhaps,  for  1  <  k  <  n  —  1  and  1  <  /  <  m  —  1, 

Pk,  o  =  Pk,  1 
Pk,m—l  =  Pk,m— 2 
PO,l  =  Pi, l 
Pn-\,l  =  Pn-2,/ 

and  similarly  for  q  (points  in  the  corner  are  copied  from  the  nearest  neighbor 
diagonally  in  the  interior  of  the  region).  It  may  be  better  to  again  use  a 
different  implementation,  where  the  values  for  points  on  the  boundary  are 
computed  from  values  at  interior  cells  that  have  the  same  “color.”  That  is, 
for  2  <  k  <  n  —  2  and  2  <  /  <  m  —  2, 

1/  X 

Pk,  0  =  ^{Pk-1,1  +Pk+l,l), 

Pk,m- 1  =  r(p*-l,m-2  +  Pfc+ 1  ,m — 2  ), 

\  (100) 
Po,i  =  +Pl,/+l)> 

1 .  x 

Pn-l,l  =  2  (Pn— 2,/— 1  "b  Pn— 2,1+1 )) 

and  similarly  for  q.  As  before,  the  comer  points,  and  one  point  on  each  side 
of  the  comer  have  to  be  copied  diagonally,  without  averaging,  since  only  one 
of  the  two  values  needed  lies  in  the  interior  of  the  region. 

6.4  Iterative  Schemes  and  Parallelism 

There  are  numerous  iterative  schemes  for  solution  of  large  sparse  sets  of  equa¬ 
tions,  among  them: 

•  Gauss-Seidel — with  replacement — sequential  update; 

•  Jacobi — without  replacement — parallel  update; 

•  Successive  Over-Relaxation; 

•  Kazmarz  relaxation; 

•  Line  relaxation. 

Successive  over-relaxation  (SOR)  makes  an  adjustment  from  the  old  value 
that  is  a  times  the  correction  computed  from  the  basic  equations.  That  is, 
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for  example, 

4?+,)  =  4?1  +  a  (4?’  -  4?’)  (101) 

where  z ^  is  the  “new”  value  calculated  by  the  ordinary  scheme  without  over¬ 
relaxation.  When  q  >  1,  this  amounts  to  moving  further  in  the  direction  of 
the  adjustment  than  suggested  by  the  basic  equations.  This  can  speed  up 
convergence,  but  also  may  lead  to  instability32.  The  Gauss- Seidel  method 
typically  can  be  sped  up  in  this  fashion  by  choosing  a  value  for  a  close  to 
two — the  scheme  becomes  unstable  for  a  >  2.  Unfortunately  the  Gauss-Seidel 
method  does  not  lend  itself  to  parallel  implementation. 

The  Jacobi  method  is  suited  for  parallel  implementation,  but  successive 
over- relaxation  cannot  be  applied  directly — the  scheme  diverges  for  a  >  1. 
This  greatly  reduces  the  speed  of  convergence.  Some  intuition  may  be  gained 
into  why  successive  over-relaxation  cannot  be  used  in  this  case,  when  it  is 
noted  that  the  neighbors  of  a  particular  cell,  the  ones  on  which  the  future  value 
of  the  cell  is  based,  are  changed  in  the  same  iterative  step  as  the  cell  itself. 
This  does  not  happen  if  we  use  the  Gauss-Seidel  method,  which  accounts  for 
its  stability.  This  also  suggests  a  modification  of  the  Jacobi  method,  where 
the  parallel  update  of  the  cells  is  divided  into  sequential  updates  of  subsets 
of  the  cells.  Imagine  coloring  the  cells  in  such  a  way  that  the  neighbors  of  a 
given  cell  used  in  computing  its  new  value  have  a  different  color  from  the  cell 
itself.  Now  it  is  “safe”  to  update  all  the  cells  of  one  color  in  parallel  (for  an 
analogous  solution  to  a  problem  in  binary  image  processing,  see  chapter  4  of 
[Horn  86]). 

Successive  over- relaxation  can  be  used  with  this  modified  Jacobi  method. 
If  local  averages  are  computed  using  only  the  four  edge- adjacent  neighbors  of 
a  cell,  then  only  two  colors  are  needed  (where  the  colors  axe  assigned  according 
to  whether  i  +  j  is  even  or  odd — see  Figure  4).  Each  step  of  the  iteration  is 
carried  out  in  two  sub-steps,  one  for  each  of  the  cells  of  one  color.  The  above 
shows  that  the  improved  convergence  rates  of  successive  over-relaxation  can 
be  made  accessible  to  parallel  implementations. 

When  the  illumination  of  the  surface  is  oblique  (light  source  away  from 
the  viewer),  R(p,q)  will  tend  to  be  locally  approximately  linear.  This  means 
that  the  gradient  of  R(p ,  q)  will  point  in  more  or  less  the  same  direction  over 
some  region  of  the  image.  The  effect  of  this  is  that  influences  on  the  adjust¬ 
ments  of  the  estimated  gradient  tend  to  be  much  smaller  along  a  direction  at 
right  angles  to  the  direction  “away  from  the  light  source,”  than  they  are  along 
other  directions.  This  can  be  seen  most  easily  when  the  coordinate  system 
is  aligned  with  the  direction  toward  a  single  light  source  in  such  a  way  that 

32Conversely,  if  the  basic  method  has  a  tendency  to  be  unstable,  then  one  can 
“under-relax” — that  is,  use  a  value  a  <  1. 
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Figure  4.  The  modified  Jacobi  method  operates  on  subsets  of  cells  with 
different  “colors”  at  different  times.  In  the  simplest  case,  there  are  only 
two  colors,  one  for  the  cells  where  the  sum  of  the  indexes  is  even,  the  other 
for  the  cells  where  the  sum  of  the  indexes  is  odd. 


the  reflectance  map  has  bilateral  symmetry  with  respect  to  the  axis  q  =  0. 
Then  Rq  will  be  small,  at  least  for  gradients  near  the  p-axis.  In  this  case  the 
coefficients  on  the  diagonal  of  the  2x2  matrix  may  be  very  different  in  mag¬ 
nitude.  This  is  analogous  to  a  system  of  equations  being  much  stiffer  in  one 
direction  than  another,  and  suggests  that  the  convergence  rate  may  be  lower 
in  this  case.  A  possible  response  to  this  difficulty  is  the  use  of  line  relaxation. 


6.5  Aliasing,  and  How  to  Avoid  It 

Discrete  samples  can  represent  a  continuous  waveform  uniquely  only  if  the  con¬ 
tinuous  waveform  does  not  contain  frequency  components  above  the  Nyquist 
rate  (wo  =  tt/c,  where  e  is  the  spacing  between  samples).  If  a  waveform  is 
sampled  that  contains  higher  frequency  components,  these  make  contributions 
to  the  sampled  result  that  are  not  distinguishable  from  low  frequency  compo¬ 
nents.  If,  for  example,  we  have  a  component  at  frequency  ujq  <  u>  <  2u>o,  it 
will  make  the  same  contributions  as  a  component  at  frequency  2u>o  —  This 
is  what  is  meant  by  aliasing.  Ideally,  the  continuous  function  to  be  sampled 
should  first  be  lowpass  filtered.  Filtering  after  sampling  can  only  suppress 
desirable  signal  components  along  with  aliased  information. 

Numerical  estimation  of  derivatives  is  weakly  ill-posed.  The  continuous 
derivative  operator  multiplies  the  amplitude  of  each  spatial  frequency  com¬ 
ponent  by  the  frequency,  thus  suppressing  low  frequencies  and  accentuating 
higher  frequencies.  Any  corruption  of  the  higher  frequencies  is  noticeable, 
particularly  if  most  of  the  signal  itself  is  concentrated  at  lower  frequencies. 
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This  means  that  we  have  to  be  careful  how  we  estimate  derivatives  and  how 
we  sample  the  image. 

Suppose,  for  example,  that  we  have  an  image  of  a  certain  size,  but  that  we 
would  like  to  run  our  shape-from- shading  algorithm  on  a  smaller  version,  per¬ 
haps  to  obtain  a  result  in  a  reasonable  amount  of  time,  or  to  cheaply  provide 
useful  initial  values  for  iteration  on  the  finer  grid.  It  would  be  quite  wrong 
to  simply  sub-sample  the  original  image.  Simple  block- averaging  is  better,  al¬ 
though  frequency  analysis  shows  that  the  response  of  a  block-averaging  filter 
first  drops  to  zero  only  at  twice  the  Nyquist  frequency.  It  is  better  to  use  a 
cubic  spline  approximation  of  the  ideal 

sin(7rx/e) 

(ttx/c) 

response  for  filtering  before  sub-sampling  [Rifman  &  McKinnon  74]  [Bern¬ 
stein  76]  [Keys  81]  [Abdou  &  Young  82].  There  is  nothing  specific  in  the  above 
relating  to  shape-from-shading;  these  are  considerations  that  apply  generally 
to  machine  vision. 

Similar  notions  apply  to  processing  of  the  surface  itself.  If  we  have  a 
digital  terrain  model  of  a  certain  resolution  and  want  to  generate  a  lower 
resolution  shaded  image  from  it,  we  need  to  first  filter  and  sample  the  digital 
terrain  model.  Otherwise  the  result  will  be  subject  to  aliasing,  and  some 
features  of  the  shaded  image  will  not  relate  in  a  recognizable  way  to  features 
of  the  surface. 

Finally,  in  creating  synthetic  data  it  is  not  advisable  to  compute  the 
surface  gradient  on  a  regular  discrete  set  of  points  and  then  use  the  reflectance 
map  to  calculate  the  expected  brightness  values.  At  the  very  least  one  should 
perform  this  computation  on  a  grid  that  is  much  finer  than  the  final  image, 
and  then  compute  block  averages  of  the  result  to  simulate  the  effect  of  finite 
sensing  element  areas — just  as  is  done  in  computer  graphics  to  reduce  aliasing 
effects33 . 

(This  hints  at  an  interesting  problem,  by  the  way,  since  the  brightness 
associated  with  the  average  surface  orientation  of  a  patch  is  typically  not  quite 
equal  to  the  average  brightness  of  the  surface,  since  the  reflectance  map  is  not 
linear  in  the  gradient.  This  means  that  one  has  to  use  a  reflectance  map 
appropriate  to  the  resolution  one  is  working  at — the  reflectance  map  depends 
on  the  optical  properties  of  the  micro-structure  of  the  surface,  and  what  is 
micro-structure  depends  on  what  scale  one  is  viewing  the  surface  at.) 

330ne  can  obtain  good  synthetic  data,  however,  with  an  exact  algebraic  solution, 
by  sampling  the  height  on  a  regular  discrete  set  of  points  and  then  estimating 
the  derivatives  numerically,  as  discussed  in  section  5.1.  This  was  done  here  to 
generate  most  of  the  examples  shown  in  section  7. 


6.  Some  Implementation  Details 
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6.6  Measuring  the  Quality  of  Reconstruction 

There  are  many  ways  of  accessing  the  quality  of  the  solution  surface  generated. 

Not  all  are  useful: 

•  In  the  case  of  a  synthetic  image  obtained  from  a  surface  model,  the  best 
test  of  the  output  of  a  shape- from-shading  algorithm  is  comparison  of  the 
surface  orientations  of  the  computed  result  with  those  of  the  underlying 
surface.  One  can  either  compute  the  root-mean-square  deviation  of  the 
directions  of  the  computed  surface  normals  from  the  true  surface  normals, 
or  just  the  root-mean-square  difference  in  the  gradients  themselves. 

•  Shading  is  a  function  of  the  surface  gradient  and  thus  most  sensitive 
to  higher  spatial  frequencies.  Conversely,  in  the  presence  of  noise  and 
reconstruction  errors,  we  expect  that  the  lower  spatial  frequencies  will  not 
be  recovered  as  well  as  the  higher  ones.  This  makes  pointwise  comparison 
of  the  heights  of  the  computed  surface  with  that  of  the  original  surface 
somewhat  less  useful,  since  errors  in  the  lower  spatial  frequencies  will 
affect  this  result  strongly.  Also,  errors  in  height  will  be  a  function  of  the 
width  of  the  region  over  which  one  has  attempted  to  recover  height  from 
gradient. 

•  Comparison  of  an  “image”  obtained  by  making  brightness  a  function 
of  height  with  a  similar  “image”  obtained  from  the  original  surface  is 
usually  also  not  very  useful,  since  such  a  representation  is  not  sensitive 
to  surface  orientation  errors,  only  gross  errors  in  surface  height.  Also, 
people  generally  find  such  displays  quite  hard  to  interpret. 

•  Oblique  views  of  “wire-meshes”  or  “block-diagrams”  defined  by  the  dis¬ 
crete  known  points  on  the  surface  may  be  helpful  to  get  a  qualitative  idea 
of  surface  shape,  but  can  be  misleading  and  are  difficult  to  compare.  If 
the  shape-from-shading  scheme  is  working  anything  like  it  is  supposed 
to,  the  differences  between  the  solution  and  the  true  surface  are  likely  to 
be  too  small  to  be  apparent  using  this  mode  of  presentation. 

•  Comparing  the  original  image  with  an  image  obtained  under  the  same 

fighting  conditions  from  the  solution  is  not  useful,  since  the  brightness 
error  is  reduced  very  quickly  with  most  iterative  schemes.  Also,  a  “so¬ 
lution”  can  have  gradient  field  that  yields  exactly  the  correct 

image  when  illuminated  appropriately,  yet  it  may  not  even  be  integrable. 
In  fact,  the  “surface”  may  yield  an  arbitrary  second  image  when  illumi¬ 
nated  from  a  different  direction  unless  p  and  q  are  forced  to  be  consistent 
(that  is,  unless  py  =  qx)  as  discussed  at  the  end  of  section  7.3. 

•  If  the  underlying  surface  is  known,  shaded  views  of  the  solution  and 
the  original  surface,  produced  under  fighting  conditions  different  from 
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those  used  to  generate  the  input  to  the  algorithm,  are  worth  comparing. 
This  is  a  useful  test,  that  immediately  shows  up  shortcomings  of  the 
solution  method.  It  also  is  a  graphic  way  of  portraying  the  progress 
of  the  iteration — one  that  is  easier  to  interpret  than  a  set  of  numbers 
representing  the  state  of  the  computation. 

•  Various  measures  of  departure  from  integrability  may  be  computed.  Per¬ 
haps  most  useful  are  comparisons  of  numerical  estimates  of  ( zx,zy )  with 
(p,g).  Slightly  iCoaj  useful  is  the  difference  ( py  —  qt)  the  solution,  since 
the  height  z  may  still  not  have  converged  to  the  best  fit  to  p  and  <7,  even 
when  the  gradient  itself  is  almost  integrable. 


6.7  When  to  Stop  Iterating 

As  is  the  case  with  many  iterative  processes,  it  is  difficult  to  decide  when  to 
stop  iterating.  If  we  knew  what  the  underlying  surface  was,  we  could  just  wait 
for  the  gradient  of  the  solution  to  approach  that  of  the  surface.  But,  other 
than  when  we  test  the  algorithm  on  synthetic  images,  we  do  not  know  what 
the  surface  is,  otherwise  we  would  probably  not  be  using  a  shape- from-shading 
method  in  the  first  place!  Some  other  tests  quantities  include: 

•  The  brightness  error 

JJ (E(x,y)  -  R(p,q)Y  dxdy  (103) 

should  be  small.  Unfortunately  this  error  becomes  small  after  just  a  few 
iterations,  so  it  does  not  yield  a  useful  stopping  criterion. 

•  The  departure  from  smoothness 

J J  {pi  +  p2y  +  ql  +  q2y)  dx  dy  (104) 

also  drops  as  the  solution  is  approached,  but  it  does  not  constitute  a  par¬ 
ticularly  good  indicator  of  approach  to  the  solution.  In  particular,  when 
one  comes  close  to  the  solution,  one  may  wish  to  reduce  the  parameters 
A,  perhaps  even  to  zero,  in  which  case  further  iterations  may  actually 
reduce  smoothness  in  order  to  better  satisfy  the  remaining  criteria. 

•  One  of  the  measures  of  lack  of  integrability 

J I {Py  ~  ?*)2  dx  dy  (105) 

is  also  not  too  useful,  since  it  can  at  times  become  small, or  stop  changing 
significantly,  even  when  z  is  still  inconsistent  with  p  and  q. 

•  Another  measure  of  lack  of  integrability 

JJ ((zx  -  p)2 +  (z9  -  q)2)  dxdy  (106) 


7.  Some  Experimental  Results 
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appears  to  be  very  useful,  since  it  drops  slowly  and  often  keeps  on  chang¬ 
ing  until  the  solution  has  converged. 

•  One  can  also  keep  track  of  the  rate  of  change  of  the  solution  with  itera¬ 
tions 

II  (i)’ *(!)’**■  ■“> 

One  should  not  stop  until  this  has  become  small.  In  most  cases  it  helps 
to  continue  for  a  while  after  the  above  measures  stop  changing  rapidly 
since  the  solution  often  continues  to  adjust. 

Some  of  the  implementation  details  given  above  may  appear  to  be  extraneous. 
However,  when  all  of  these  matters  are  attended  to,  then  the  iterative  algo¬ 
rithm  will  not  “walk  away”  from  the  solution,  and  it  will  find  the  solution, 
to  machine  precision,  given  exact  data  (and  assuming  that  boundary  condi¬ 
tions  for  p  and  q  are  given,  and  that  A'  is  reduced  to  zero  as  the  solution  is 
approached).  Convergence  to  the  exact  solution  will  not  occur  when  some¬ 
thing  is  amiss,  such  as  a  mismatch  between  the  discrete  estimators  of  the  first 
derivative  and  the  discrete  estimator  of  the  Laplacian.  It  is  not  yet  clear  how 
significant  all  of  this  is  when  one  works  with  real  image  data,  where  there  is  no 
exact  solution,  and  where  the  error  introduced  by  incorrect  implementation 
detail  may  be  swamped  by  errors  from  other  sources. 

7.  Some  Experimental  Results 

The  new  algorithm  has  been  applied  to  a  number  of  synthetic  images  of  simple 
shapes  (such  as  an  asymmetrical  Gaussian,  a  sum  of  Gaussian  blobs,  and  a 
sum  of  low  frequency  sinusoidal  gratings)  generated  with  a  number  of  different 
reflectance  maps  (including  one  linear  in  p  and  q,  Lambertian  with  oblique 
illumination,  and  a  rotationally  symmetric  one).  These  synthetic  images  were 
small  (usually  64  x  64  picture  cells)  in  order  to  keep  the  computational  time 
down.  Typically  the  surface  normals  would  be  within  a  degree  or  two  of  the 
correct  direction  after  a  few  hundred  iterations.  With  appropriate  boundary 
conditions,  the  computed  shape  would  eventually  be  the  same  (to  machine 
precision)  as  the  shape  used  to  generate  the  synthetic  image.  In  each  case, 
the  brightness  error  decreased  rapidly,  while  the  integrability  of  the  estimated 
gradient  decreased  much  more  slowly. 


7.1  Graphical  Depiction  of  Solution  Process 

For  help  in  debugging  the  algorithm,  and  for  purposes  of  determining  a  good 
schedule  for  adjusting  the  parameters  p  and  A',  it  is  useful  to  print  out  the 
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(run-schedule  2) 

Lanbda:  l.eBBB  flu:  8.1888  Eps:  1.8888  Rlpha-pq;  1.7808  Rlpha-z:  1.7888  Cor-pq  8.2588  Cor-t  1.8888  (64  X  64) 

Iter:  8  Grad  E:  2.09619  Bright  E:  8.52570  Int-z  E:  0.00000  (Int-pq  E:  8.00000)  Unsnooth:  3.24605 

Iter:  4  Grad  E:  0.60296  Bright  E:  0.17914  Int-z  E:  1.38119  (Int-pq  E:  0.47459)  Unsnooth:  0.66932  dpq:l. 091998  dt:2. 173221 

Iter:  8  Grad  E:  8.27183  8right  E:  0.07793  Int-z  E:  0.32173  (Int-pq  E:  8.89847)  Unsnooth:  0.19669  dpq:0. 292286  dz:0.522?6S 

Lanbda:  3.5008  Hu:  0.1008  Eps:  1.0000  81pha-pq:  1.7888  Rlpha-z:  1.7880  Cor-pq  0.2588  Cor-z  1.8880  (64  X  64) 

Iter:  14  Grad  E:  8.21348  Bright  E:  0.03485  Int-z  E:  0.88287  (Int-pq  E:  8.82512)  Unsnooth:  0.06240  doq:0. 044399  dz:0. 870357 

Iter:  20  Grad  E:  0.20982  Bright  E:  0.03317  Int-z  E:  8.86199  (Int-pq  E:  0.31987)  Unsnooth:  0.05529  dpq:0. 006952  dr: 0.02401 9 

Iter:  26  Grad  E:  0.28821  Bright  E:  0.83297  Int-z  E:  0.05496  (Int-pq  El  0.01918)  Unsnooth:  0.0S508  dpq:0. 002019  dz:0. 016154 

Iter:  32  Grad  E:  0.20702  Bright  E:  0.83268  Int-z  E:  0.05077  (Int-pq  E:  0.01886)  Unsnooth:  0.05510  dpq:0. 001203  dz:8. 012275 

Lanbda:  0.2000  Mu:  0.1000  Eps:  1.0000  Rlpha-pq:  1.7880  Rlpha-z;  1.7808  Cor-pq  8.2580  Cor-z  1.0088  (64  X  64) 

Iter:  40  Grad  E:  0.187B5  Bright  E:  0.82212  Int-z  E:  0.04226  (Int-pq  E:  0.02146)  Unsnooth:  0.06969  dpq:8. 082513  dz:0. 011044 

Iter:  48  Grad  E:  0.18496  Bright  E:  0.02183  Int-z  E:  0.03885  (Int-pq  E:  0.02111)  Unsnooth:  0.07009  dpq:0. 000822  dz:0. 007575 

Iter:  56  Grad  E:  0.18334  Bright  E:  0.02174  Int-z  £:  0.03689  (Int-pq  E:  0.02091)  Unsnooth:  0.07019  dpq:0. 880601  dz:0. 006001 

Iter:  64  Grad  E:  0.18226  Bright  E:  0.02169  Int-z  E:  0.03562  (Int-pq  E:  0.02077)  Unsnooth:  0.07022  dpq:0. 000485  dz:0. 005806 

Lanbda:  0.1000  hu:  0.1000  Eps:  1.0000  Rlpha-pq:  1.7000  Rlpha-z:  1.7000  Cor-pq  0.2500  Cor-z  1.0080  (64  X  64) 

Iter:  88  Grad  E:  0.16684  Bright  E:  0.01523  Int-z  E:  0.02934  (Int-pq  E:  0.02160)  Unsnooth:  0.08152  dpq:0. 008535  dz:0. 004017 

Iter:  96  Grad  E:  0.16446  Bright  E:  0.01516  Int-z  E:  0.02B01  (Int-pq  E:  0.02143)  Unsnooth:  0.08167  dpq:0. 000324  dz:0. 082881 

Iter:  112  Grad  E:  0.16381  Bright  £:  0.01514  Int-z  E:  0.02736  (Int-pq  E:  0.02136)  Unsnooth:  0.08170  dpq:0. 808252  dz:0. 002305 

Iter:  128  Grad  E:  0.16195  Bright  E:  3.01513  Int-z  E:  0.02696  vlnt-pq  E:  3.02131)  Unsnooth:  0.08170  dpq:0. 008204  dz:0. 801920 

Lanbda:  8.0500  hu:  0.1800  Eps:  1.0000  Rlpha-pq:  1.7000  Rlpha-z:  1.7000  Cor — pq  3.2500  Cor — z  1.0608  (64  X  64) 

Iter:  168  Grad  £:  0.14655  Bright  £:  0.81822  Int-z  E:  0.02194  (Int-pq  E:  0.02062)  Unsnooth:  0.09233  doq:0. 000219  dz:0. 001501 

Iter:  192  Grad  E:  0.14411  Bright  E:  0.01821  Int-z  E:  0.02155  (Int-pq  E:  0.02056)  Unsnooth:  0.09246  doq:0. 000142  dz:0. 001074 

Iter:  224  Grad  E:  0.14252  Bright  £:  0.01021  Int-z  E:  8.02138  (Int-pq  £:  8.02852)  Unsnooth:  0.09253  dpq:0. 800107  dz:0. 000839 

Iter:  256  Grad  E:  0.14139  Bright  E:  0.01021  Int-z  E:  0.0212B  (Int-pq  E:  8.02049)  Unsnooth:  0.09257  dpq:0. 000084  dz:0. 800673 

Lanbda:  0.0200  hu:  0.0808  Eps:  1.8880  Rlpha-pq:  1.7000  Rlpha-z:  1.7000  Cor-pq  0.2508  Cor-z  1.0000  (64  X  64) 

Iter:  320  Grad  E:  0.12432  Bright  E:  0.00573  Int-z  E:  0.01760  (Int-pq  E:  0.01893)  Unsnooth:  8.10434  dpq:0. 000101  dz:0. 880698 

Iter:  384  Grad  E:  0.11874  8right  E:  0.00394  Int-z  E:  0.01262  (Int-pq  E:  0.01359)  Unsnooth:  0.11303  dpq:0. 000586  dz:0. 001413 

Iter:  448  Grad  E:  0.10616  Bright  E:  8.00426  Int-z  E:  0.01341  (Int-pq  E:  0.01433)  Unsnooth:  0.11413  dpq:0. 810938  dr: 0.004026 

Iter:  512  Grad  E:  0.18092  Bright  E:  0.00381  Int-z  E:  8.01180  (Int-pq  E:  0.01326)  Unsnooth:  0.11442  dpq:0. 000078  dr: 0.000526 

Lanbda:  0.0100  hu:  0.0500  Eps:  1.0080  Rlpha-pq:  1.7000  Rlpha-z:.  1.7000  Cor-pq  0.2500  Cor-z  1.8800  (64  X  64) 

Iter:  640  Grad  E:  0.09924  Bright  E:  0.00378  Int-z  E:  8.01167  (Int-pq  £:  0.81318)  Unsnooth:  0.11467  dpq:0. 000068  dr: 0.000432 

Iter:  768  Grad  E:  0.09640  Bright  E:  0.00376  Int-z  E:  0.01153  (Int-pq  E:  0.01309)  Unsnooth:  0.11514  dpq:8. 000033  dr: 0.000251 

Iter:  896  Grad  E:  0.09519  Bright  E:  0.00374  Int-z  E:  0.01147  (Int-pq  E:  0.81302)  Unsnooth:  0.11537  dpq:0. 000021  dz:0.00«>164 

Iter:1024  Grad  E:  0.09459  Bright  E:  8.08373  Int-z  E:  0.01144  (Int-pq  E:  8.01299)  Unsnooth:  8.11551  dpq:0. 800014  dz:8. 000111 
Lanbda:  0.0050  hu:  8.0380  Eps:  1.0000  Rlpha-pq:  1.7800  Rlpha-z:  1.7008  Cor-pq  0.2500  Cor-z  1.8800  (64  X  64) 

Iter:1280  Grad  E:  0.10298  Bright  E:  8.00202  Int-z  E:  0.01598  (Int-pq  E:  0.01850)  Unsnooth:  0.11612  dpq : 0 . 003722  dz:0. 000633 

I ter: 1536  Grad  E:  0.10330  Bright  E:  0.00201  Int-z  E:  0.01592  (Int-pq  E:  0.01841)  Unsnooth:  8.11608  dpq:0. 003155  dz:0. 000624 

I tar: 1792  Grad  E:  0.10336  Bright  E:  8.00201  Int-z  E:  0.01590  (Int-pq  E:  0.01839)  Unsnooth:  8.11608  dpq:0. 003207  dr: 0.008635 

Iter : 2048  Grad  E:  0.10340  Bright  E;  0.00201  Int-z  E:  0.01590  (Int-pq  E:  0.01839)  Unsnooth:  0.11607  dpq:0. 803141  dz: 0.800629 

Lanbda:  0.0020  Mu:  0.0200  Eps:  1.0000  Rlpha-pq:  1.7000  Rlpha-z:  1.7800  Coi — pq  0.2500  Cor — z  1.0000  (64  X  64) 

I  ter : 2304  Grad  E:  0.09030  Bright  E:  0.00111  Int-z  E:  0.01341  (Int-pq  E:  0.01576)  Unsnooth:  8.12242  dpq:0. 007526  dz:0. 001322 

Iter: 2560  Grad  E:  0.08954  Bright  E:  0.00110  Int-z  E:  0.01337  (Int-pq  E:  0.01570)  Unsnooth:  8.12259  dpq:0. 007359  dz:0. 001289 

Iter:2016  Grad  E:  0.08932  Bright  E:  8.80110  Int-z  E:  0.01336  (Int-pq  E:  0.01568)  Unsnooth:  8.12264  dpq:0. 007298  dz:0. 001277 

Iter: 3072  Grad  E:  0.08924  Bright  E:  0.00110  Int-z  E:  0.01336  (Int-pq  E:  0.01567)  Unsnooth:  8.12267  dpq:0. 007276  dz:0. 001273 

Lanbda:  8.0010  Mu:  0.0100  Eps:  1.0000  Rlpha-pq:  1.7000  Rlpha-z:  1.7000  Cor-pq  0.2580  Cor-z  1.0000  (64  X  64) 

I ter: 3328  Grad  E:  0.08869  Bright  E:  8.00086  Int-z  E:  0.01434  (Int-pq  E:  0.01679)  Unsnooth:  8.12387  dpq:0. 011513  dz:0. 002398 

I ter: 3584  Grad  E:  0.08866  Bright  E:  0.00086  Int-Z  E:  0.01433  (Int-pq  E:  0.01679)  Unsnooth:  0.12388  dpq:0. 011497  dz:0. 002395 

I ter: 3840  Grad  E:  0.08865  Bright  E:  0.80086  Int-z  E:  0.01433  (Int-pq  E:  0.01678)  Unsnooth:  0.12388  dpq:0. 811493  dz:0. 002394 

I ter: 4096  Grad  E:  0.08864  Bright  E:  0.00086  Int-z  E:  0.01433  (Int-pq  E:  0.01678)  Unsnooth:  0.12388  dpq:0. 811492  dz:0. 002394 

Lanbda:  0.0005  Mu:  0.0080  Eps:  1.0000  Rlpha-pq:  1.7000  Rlpha-z:  1.7000  Co: — pq  0.2500  Co: — z  1.0000  (64  X  64) 

I  ter : 4352  Grad  E:  8.07698  Bright  E:  0.00050  Int-z  E:  0.01111  (Int-pq  E:  0.01290)  Unsnooth:  0.12855  dpq:0. 007626  dz:0. 001396 

Iter: 4608  Grad  E:  0.07620  Bright  E:  0.00049  Int-z  E:  0.01108  (Int-pq  E:  0.01285)  Unsnooth:  0.12875  dpo:0. 007546  dz:0. 001378 

I ter: 4864  Grad  E:  0.07597  Bright  E:  0.00049  Int-z  E:  0.01187  (Int-pq  E:  0.01284)  Unsnooth:  0.12882  dpq:0. 807509  dz:0. 001371 

I ter: 51 20  Grad  E:  0.07588  Bright  E:  0.00049  Int-z  E:  0.01107  (Int-pq  E:  0.01283)  Unsnooth:  0.12885  dpq:0. 807494  dz:0. 001368 

Lanbda:  0.0003  Mu:  0.0070  Eps:  1.0000  Rlpha-pq:  1.7000  Rlpha-z:  1.7000  Co: — pq  0.2500  Co: — z  1.0000  (64  X  64) 

I  ter: 5376  Grad  E:  0.06687  Bright  E:  0.00030  Int-z  E:  0.00890  (Int-pq  E:  0.01026)  Unsnooth:  0.13258  dpq:0. 005249  dz:0. 000933 

Iter: 5632  Grad  E:  0.06620  Bright  E:  0.00030  Int-z  E:  0.00887  (Int-pq  E:  0.01022)  Unsnooth:  0.13278  dpq:0. 805163  dz:0. 008917 

I  ter: 5088  Grad  £:  0.06599  Bright  E:  0.00030  Int-z  E:  0.00887  (Int-pq  E:  8.01021)  Unsnooth:  8.13285  dpq:0. 005119  dz:0. 008909 

I ter: 61 44  Grad  E:  0.06590  Bright  E:  8.00030  Int-z  E:  0.00886  (Int-pq  E:  0.01020)  Unsnooth:  0.13288  dpq : 0 . 805099  dz:0. 000986 

Lanbda:  0.0002  Mu:  0.0058  Eps:  1.0000  Rlpha-pq:  1.7000  Rlpha-z:  1.7000  Cor-pq  8.2508  Cor-z  1.0000  (64  X  64) 

I  ter : 6656  Grad  E:  0.06406  Bright  E:  0.00025  Int-z  E:  8.08861  (Int-pq  E:  0.08988)  Unsnooth:  0.13378  dpq:0. 005282  dz: 0.000968 

Iter: 7168  Grad  E:  0.06399  Bright  E:  8.00025  Int-z  E:  0.00860  (Int-pq  E:  0.00988)  Unsnooth:  0.13381  dpq:0. 005270  dz:0. 008966 

I  ter : 7680  Grad  E:  0.06398  Bright  E:  0.80025  Int-z  E:  0.00860  (Int-pq  E:  0.00988)  Unsnooth:  0.13381  dpq:0. 085268  dz:0. 008965 

I ter: 81 92  Grad  E:  0.06397  Bright  E:  0.80025  Int-z  E:  0.00860  (Int-pq  E:  0.00988)  Unsnooth:  0.13382  dpq:0. 005268  dz:0. 000965 

Lanbda:  8.0001  Mu:  0.0020  Eps:  1.0000  Rlpha-pq:  1.7000  Rlpha-z:  1.7000  Cor-pq  0.2500  Cor-z  1.0000  (64  X  64) 

Iter:8704  Grad  E:  8.06941  Bright  E:  0.00032  Int-z  E;  0.01012  (Int-pq  E:  0.01162)  Unsnooth:  0.13186  dpq:0. 007721  dz:0. 001530 

Iter:9216  Grad  E:  0.06955  Bright  E:  0.00032  Int-z  E:  0.01012  (Int-pq  E:  0.  J1163)  Unsnooth:  0.13182  dpq:0. 007740  dz:0. 001535 

I  ter: 9728  Grad  E:  0.06958  Bright  E:  8.00032  Int-z  E:  0.01012  (Int-pq  E:  0.01163)  Unsnooth:  0.13181  dpq:0. 007744  dz:0. 001536 

Iter : 10240  Grad  E:  0.069S8  Bright  E:  0.80032  Int-z  E:  0.81012  (Int-pq  E:  0.01163)  Unsnooth:  0.13180  dpq:0. 007744  dz:8. 001536 
Lanbda:  0.0000  Mu:  0.0010  Eps:  1.0000  R)pha-pq:  1.7000  Rlpha-z:  1.7000  Cor-pq  8.2500  Cor-z  1.0000  (64  X  64) 

Iter:10752  Grad  E:  0.01967  Bright  E:  0.80081  Int-z  E:  0.80106  (Int-pq  E:  0.00092)  Unsnooth:  0.16893  dpq:0. 000062  dz:0. 000185 

Iter : 11264  Grad  E:  0.00742  Bright  E:  0.00000  Int-z  E:  0.00039  (Int-pq  E:  0.00038)  Unsnooth:  0.16351  dpq:0. 000013  dz:8. 000062 

Iter:11776  Grad  E:  0.00366  Bright  E:  0.00000  Int-z  E:  0.00017  (Int-pq  E:  0.00014)  Unsnooth:  0.16446  dpq:0. 000085  dz:0. 000025 

I  ter : 12288  Grad  E:  0.00189  Bright  E:  0.00000  Int-z  E:  0.08008  (Int-pq  E:  0.00007)  Unsnooth:  0.16488  dpq:8. 000803  dz: 0.000011 


Figure  5.  Diagnostic  trace  of  various  error  measures.  This  sequence  of 
results  corresponds  to  the  reconstruction  of  the  sharp-edged  crater  shape 
shown  in  Figure  7.  This  kind  of  presentation  is  important,  but  must  be 
supplemented  by  some  graphic  depiction  of  the  evolving  solution  surface. 
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diagnostic  measurements  discussed  in  sections  6.6  and  6.7.  But  it  is  hard 
to  tell  exactly  what  is  going  on  just  by  looking  at  a  large  table  of  numbers 
such  as  that  shown  in  Figure  5.  It  is  important  to  also  provide  some  graphic 
depiction  of  the  evolving  shape  as  it  is  computed.  To  make  shaded  images 
of  the  reconstructed  surface  useful,  however,  they  must  be  illuminated  from 
a  direction  different  from  the  direction  of  illumination  used  for  the  original 
input  image34.  Shown  in  Figure  6,  is  such  a  sequence  of  shaded  images  gener¬ 
ated  during  the  reconstruction  of  the  surface  of  a  polyhedral  object,  starting 
from  a  random  field  of  surface  orientations.  Here  the  image  provided  to  the 
algorithm35  corresponded  to  illumination  form  the  Northwest,  while  illumi¬ 
nation  from  the  Northeast  was  used  to  display  the  reconstruction.  Note  how 
the  edges  become  sharper  as  A',  controlling  the  contribution  of  the  penalty 
term  for  departure  from  smoothness,  is  made  smaller  and  smaller.  This  ex¬ 
ample  also  illustrates  the  algorithm’s  ability  to  deal  with  surfaces  that  have 
discontinuities  in  surface  orientation. 

Because  of  the  interest  in  application  to  astrogeology,  a  crater-like  shape 
was  also  reconstructed,  as  shown  in  Figure  7.  In  this  case,  the  algorithm 
rapidly  found  a  shape  that  was  generally  correct,  except  for  flaws  in  places  on 
the  rim  of  the  crater  in  the  Northwest  and  Northeast.  These  are  areas  where 
there  is  little  contrast  between  the  inside  and  the  outside  of  the  crater  in  the 
input  image36.  It  took  the  algorithm  a  considerable  number  of  additional 
iterations  to  determine  the  correct  continuation  of  the  shape  computed  in 
other  image  areas. 

7.2  Emergent  Global  Organization 

Often  progress  towards  the  correct  solution  is  not  as  uneventful.  Frequently 
small  internally  consistent  solution  patches  will  establish  themselves,  with 
discontinuities  in  surface  orientation  where  these  patches  adjoin.  Also,  coni¬ 
cal  singularities  form  that  tend  to  move  along  the  boundaries  between  such 
regions  as  the  iterative  solution  progresses.  Conversely,  boundaries  between 
solution  patches  often  form  along  curves  connecting  conical  singularities  that 
form  earlier.  After  a  large  enough  number  of  iterations,  patches  of  local  orga¬ 
nization  tend  to  coalesce  and  lead  to  emergent  global  organization.  This  can 

34The  test  illumination  should  be  quite  different  from  the  illumination  used  to 
generate  the  original  image — preferrably  lying  in  a  direction  that  differs  from  the 
original  source  direction  by  as  much  as  tt/2 
35The  input  image  is  not  shown,  but  is  just  like  the  last  image  in  the  sequence 
shown,  except  that  left  and  right  are  reversed. 

36 Again,  the  input  image  is  not  shown,  but  is  like  the  last  image  in  the  sequence 
shown,  except  that  left  and  right  are  reversed. 
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Figure  7.  Reconstruction  of  a  crater-like  shape.  Points  on  the  rim  in 
the  Northeast  and  the  Southwest  correspond  to  places  in  the  input  image 
where  there  is  least  contrast  between  the  inside  and  the  outside,  since  the 
direction  of  the  incident  illumination  is  parallel  to  the  rim  there. 
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be  observed  best  when  A'  is  smaller  than  what  it  would  normally  be  set  to  for 
rapid  convergence.  In  Figures  8  and  9,  for  example,  are  shown  a  sequence  of 
shapes  leading  finally  to  a  spherical  cap  on  a  planar  surface.  Within  some  re¬ 
gions,  solution  surface  patches  quickly  establish  themselves  that  individually 
provide  good  matches  to  corresponding  parts  of  the  input  image.  The  bor¬ 
ders  between  these  internally  consistent  regions  provide  error  contributions 
that  the  algorithm  slowly  reduces  by  moving  the  boundaries  and  incremen¬ 
tally  changing  the  shapes  within  each  of  the  regions.  Too  rapid  a  reduction 
of  A'  can  remove  the  incentive  to  reduce  the  creases  and  kinks  and  to  freeze 
the  solution  in  a  state  where  some  unnecessary  discontinuities  remain.  If,  for 
example,  A'  were  to  be  set  to  zero  with  a  “solution”  consisting  of  a  spherical 
cap  with  an  inner  disk  inverted,  as  in  the  right  hand  image  of  the  middle 
row  of  Figure  9,  there  would  be  no  incentive  to  further  reduce  the  length  of 
the  circular  discontinuity,  and  the  smooth  solution  for  this  part  of  the  image 
would  not  be  found. 

The  algorithm  was  also  applied  to  impossible  shaded  images.  Suppose, 
for  example,  that  we  are  dealing  with  a  Lambertian  surface  illuminated  by 
a  source  near  the  viewer  and  that  there  is  a  dark  smudge  in  the  middle  of  a 
large  planar  region  facing  us  (which  appears  brightly  lit).  It  turns  out  that 
there  is  no  surface  with  continuous  first  derivatives  that  could  give  rise  to  a 
shaded  image  with  a  simply  connected,  bounded  dark  region  in  the  middle 
of  a  bright  region  [Szeliski  &  Horn  89].  In  Figure  10  we  see  what  happens 
when  the  algorithm  attempts  to  find  a  solution.  Patches  grow  within  which 
the  solution  is  consistent  with  the  image,  but  which  yield  discontinuities  at 
boundaries  between  patches.  Conical  singularities  sit  astride  these  boundaries. 
For  all  random  initial  conditions  tried,  the  algorithm  eventually  eliminates  all 
but  one  of  these  conical  singularities.  The  computed  surface  is  in  fact  a 
“solution,”  if  one  is  willing  to  allow  such  singularities. 

The  graphical  method  of  presenting  the  progress  of  the  iterative  solutions 
illustrated  above  was  very  helpful  in  debugging  the  program  and  in  determin¬ 
ing  reasonable  schedules  for  reduction  of  the  parameters  A'  and  ju.  Shown  in 
Figure  11  are  some  examples  of  what  happens  when  things  go  wrong.  In  the 
top  row  are  shown  instabilities  arising  in  the  solution  for  the  crater-like  shape, 
near  the  points  where  there  is  low  contrast  between  the  inside  and  the  outside 
of  the  crater — that  is,  where  there  is  no  local  evidence  for  curvature.  These 
instabilities  can  be  suppressed  by  reducing  A'  more  slowly.  In  the  middle  row 
are  shown  patterns  resulting  from  various  programming  errors.  Finally,  in  the 
bottom  row  is  shown  the  propagation  of  an  instability  from  a  free  boundary 
when  A'  is  set  to  zero.  It  appears  that  the  process  is  not  stable  when  there 
are  neither  boundary  conditions  for  the  height  nor  for  the  gradient. 
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Figure  8.  Emergent  global  organization  of  local  nonlinear  iterative  pro¬ 
cess.  Internally  consistent  solutions  arise  in  certain  image  patches  with 
discontinuities  at  the  borders  between  regions.  The  boundaries  between 
these  patches  move,  and  the  solutions  within  the  patches  adjust  in  order 
to  reduce  the  sum  of  the  error  terms. 
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Figure  9.  Emergent  global  organization  of  local  nonlinear  iterative  pro¬ 
cess.  Neighboring  patches  coalesce,  as  conical  singularities  are  absorbed 
by  coalescing  with  other  singularities,  or  by  being  pushed  towards  a  con¬ 
tour  where  the  surface  orientation  is  discontinuous.  The  final  solution  is 


a  spherical  cap  resting  on  a  plane. 
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Figure  10.  What  happens  when  the  algorithm  is  confronted  with  an 
impossible  shaded  image?  The  input  image  here  (not  shown)  is  a  circularly 
symmetric  dark  smudge  in  a  uniformly  bright  surround.  The  light  source 
is  assumed  to  be  near  the  viewer.  The  algorithm  finds  a  “solution”  with 
a  single  conical  singularity. 
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Figure  11.  Graphical  depiction  of  instabilities  and  the  effects  of  pro¬ 
gramming  errors.  In  the  top  row  are  shown  instabilities  resulting  from  too 
rapid  reduction  of  the  penalty  term  for  departure  from  smoothness.  The 
middle  row  shows  the  results  of  various  programming  errors.  The  bottom 
row  shows  waves  of  instability  propagating  inwards  from  a  free  boundary. 
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Figure  12.  New  one-step  shape-from-shading  algorithm.  See  text! 
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In  the  past,  shape-from-shading  algorithms  have  often  been  “tested”  by 
verifying  that  the  computed  gradient  field  actually  generates  something  close 
to  the  given  input  image.  To  show  just  how  dangerous  this  is,  consider  Fig¬ 
ure  12,  which  demonstrates  a  new  non-iterative  method  for  recovering  a  “sur¬ 
face”  given  a  shaded  image.  In  Figure  12(a),  is  the  input  to  the  algorithm, 
while  Figure  12(c)  is  what  the  gradient  field  that  is  constructed  by  this  al¬ 
gorithm  looks  like  when  illuminated  in  the  same  way  as  the  original  surface. 
Now  Figure  12(b)  shows  what  the  original  surface  looks  like  when  illuminated 
from  another  direction.  As  a  test,  we  should  check  that  the  computed  gradi¬ 
ent  field  looks  the  same  under  these  illuminating  conditions.  But  behold,  it 
does  not!  In  Figure  12(d)  we  see  what  we  get  when  we  use  the  other  illumi¬ 
nating  condition.  The  “trick”  here  is  that  the  problem  of  shape  from  shading 
is  heavily  underconstrained  if  we  are  only  computing  a  gradient  field  and  not 
enforcing  integrability.  There  are  many  solutions  and  we  can,  in  fact,  impose 
additional  constraints.  The  underlying  gradient  field  here  was  computed  by 
solving  the  photometric  stereo  equations  [Woodham  78,  79,  80,  89]  for  the  two 
images  in  Figures  12(c)  and  (d)  under  the  two  assumed  lighting  conditions. 

The  new  algorithm  has  also  been  applied  to  synthetic  images  generated 
from  more  complicated  surfaces  such  as  digital  terrain  models  (DTMs)  used 
earlier  in  research  on  interpretation  of  satellite  images  of  hilly  terrain  [Horn  & 
Bachman  78]  [Sjoberg  &:  Horn  83]  and  in  automatic  generation  of  shaded  over¬ 
lays  for  topographic  maps  [Horn  81].  These  synthetic  images  were  somewhat 
larger  (the  one  used  for  Figure  1  is  of  size  231  x  178,  for  example).  In  this 
case,  the  simple  algorithm,  presented  in  section  5.3,  using  a  regularizing  term 
would  often  get  trapped  in  a  local  minimum  of  the  error  function  after  a  small 
number  of  iterations,  while  the  modified  algorithm  presented  in  section  5.6, 
exploiting  the  linearization  of  the  reflectance  map,  was  able  to  proceed  to  a 
solution  to  machine  precision  after  a  few  thousand  iterations.  Most  of  the 
surface  normals  typically  were  already  within  a  degree  or  so  of  the  correct 
direction  after  a  few  hundred  iterations. 

The  closeness  of  approach  to  the  true  solution  depends  on  several  of  the 
implementation  details  discussed  earlier.  In  particular,  it  was  helpful  to  use 
the  old  values  of  p  and  q  for  the  reference  point  in  the  linearization  of  R(p,  q ), 
rather  than  any  of  the  other  choices  suggested  earlier.  Also,  it  helps  to  use 
the  diagonal  averaging  scheme  in  the  iteration  for  height  rather  than  the  one 
based  on  edge- adjacent  neighbors. 


7.3  Rating  the  Difficulty  of  Shape-from-Shading  Problems 

Experiments  with  synthetic  shaded  images  suggests  that  certain  shape-from- 
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shading  problems  are  relatively  easy,  while  others  are  quite  difficult.  First 
of  all,  ba33o-relievoz 7  surfaces  (those  with  only  low  slopes)  are  easy  to  deal 
with  (see  also  section  2.6)  in  comparison  with  alto-relievo  surfaces  (those  with 
steep  slopes).  The  digital  terrain  model  used  for  the  experiment  illustrated 
in  Figure  1  falls  in  the  latter  category,  since  the  sides  of  the  glacial  cirque  are 
steep  and  the  individual  gullies  steeper  still. 

Typically  the  brightness  of  a  surface  patch  increases  the  more  it  is  turned 
towards  the  light  source.  If  it  is  turned  too  far,  however,  it  becomes  so  steep 
that  its  brightness  once  again  decreases.  There  is  a  qualitative  difference 
between  shape-from-shading  problems  where  none  of  the  surface  patches  are 
turned  that  far,  and  those  where  some  surface  patches  are  so  steep  as  to  have 
reduced  brightness.  In  the  latter  case,  there  appears  to  be  a  sort  of  two- 
way  ambiguity  locally  about  whether  a  patch  is  dark  because  it  has  not  been 
turned  enough  to  face  the  light  source  or  whether  it  has  been  turned  too  far. 
This  ensures  that  simplistic  schemes  will  get  trapped  in  local  minima  where 
patches  of  the  solution  have  quite  the  wrong  orientation.  Similarly,  the  more 
sophisticated  scheme  described  here  takes  many  more  iterations  to  unkink  the 
resulting  creases. 

The  transition  between  the  two  situations  depends  on  where  the  light 
source  is.  The  difficulty  is  reduced  when  the  illumination  is  oblique  (see  also 
section  2.6).  Conversely,  the  problem  is  more  severe  when  the  light  source  is 
at  the  viewer,  in  which  case  brightness  decreases  with  slope  independent  of 
the  direction  of  the  surface  gradient.  This  explains  why  the  algorithm  took 
longer  to  find  the  solution  in  the  case  of  the  spherical  cap  (Figure  8)  since 
it  was  illuminated  by  a  source  near  the  viewer.  Tt  was  more  straightforward 
to  find  the  solutions  for  the  truncated  hexahedron  and  the  crater-like  surface 
(Figures  6  and  7),  both  of  which  where  illuminated  obliquely.  The  above 
dichotomy  is  related  to  another  factor:  problems  where  the  relevant  part  of 
the  reflectance  map  is  nearly  linear  in  gradient  are  considerably  easier  to  deal 
with  than  those  in  which  the  reflectance  map  displays  strong  curvatures  of 
iso-brightness  contours. 

Smooth  surfaces,  particularly  when  convex,  can  be  recovered  easily.  Sur¬ 
faces  with  rapid  undulations  and  wrinkles,  such  as  the  digital  terrain  model 
surface  (Figure  1)  axe  harder.  Discontinuities  in  surface  orientation  are  even 
more  difficult  to  deal  with.  Note  that,  with  the  exception  of  the  digital  terrain 
model,  all  of  the  examples  given  here  involve  surfaces  that  have  some  curves 
along  which  the  surface  orientation  is  not  continuous.  The  spherical  cap,  for 
example,  lies  on  a  planar  surface,  with  a  discontinuity  in  surface  orientation 

37For  more  regarding  the  terms  basso-relievo,  mezzo-relievo  and  alto-relievo ,  see 
(Koenderink  &  van  Doom  80}. 
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where  it  touches  the  plane. 

Problems  where  boundary  conditions  are  not  available,  and  where  there 
are  no  occluding  boundaries  or  singular  points,  are  ill-posed.  Not  too  surpris¬ 
ingly  these  tend  to  lead  to  instabilities  in  the  algorithm,  particularly  when 
one  attempts  to  reduce  the  penalty  term  for  departure  from  smoothness.  In 
these  cases  instabilities  can  be  damped  out  to  some  extent  by  enforcing  the 
image  irradiance  equation  on  the  boundary  by  iterative  adjustment  of  the 
gradient  computed  from  the  discrete  approximation  of  the  natural  boundary 
conditions  for  p  and  q.  But  results  have  not  been  promising  enough  to  discuss 
here  in  more  detail. 

The  number  of  iterations  to  converge  to  a  good  solution  appears  to  grow 
almost  quadratically  with  image  size  (number  of  rows  or  columns).  This  is 
because  some  effects  have  to  “diffuse”  across  the  image.  This  means  that  the 
total  amount  of  computation  grows  almost  with  the  fourth  power  of  the  (lin¬ 
ear)  image  size.  It  is  well  known  that  ordinary  iterative  schemes  for  solving 
elliptic  partial  differential  equations  quickly  damp  out  higher  spatial  frequency 
errors,  while  low  frequency  components  are  removed  very  slowly.  The  way  to 
deal  with  this  problem  is  to  use  computation  on  coarser  grids  to  reduce  the 
low  spatial  frequency  components  of  the  error.  This  is  the  classic  multigrid 
approach  [Brandt  77,  80,  84].  It  is  clear  that  a  true  multigrid  implementation 
(as  opposed  to  a  simple  pyramid  scheme)38  would  be  required  to  pursue  this 
approach  further  on  larger  images.  This  is  mostly  to  cut  down  on  the  com¬ 
putational  effort,  but  can  also  be  expected  to  reduce  even  further  the  chance 
of  getting  caught  in  a  local  minimum  of  the  error  function.  Implementation, 
however,  is  not  trivial,  since  the  equations  are  nonlinear,  and  because  there 
are  boundary  conditions.  Both  of  these  factors  complicate  matters,  and  it  is 
known  that  poor  implementation  can  greatly  reduce  the  favorable  convergence 
rate  of  the  basic  multigrid  scheme  [Brandt  77,  80,  84]. 

Alternatively,  one  can  apply  so-called  direct  methods  for  solving  Poisson’s 
equations  [Simchony,  Chellappa  &  Shao  89]. 

8.  Conclusion 

The  original  approach  to  the  general  shape- from- shading  problem  requires 
numerical  solution  of  the  characteristic  strip  equations  that  arise  from  the 
first-order  nonlinear  partial  differential  equation  that  relates  image  irradiance 
to  scene  radiance  [Horn  70,  75].  Variational  approaches  to  the  problem  instead 

38A  naive  approach  has  one  solve  the  equations  on  a  coarse  grid  first,  with  the 
results  used  as  initial  conditions  for  a  finer  grid  solution  after  interpolation.  True 
multigrid  methods  are  more  complex,  but  also  have  much  better  properties. 
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minimize  the  sum  of  the  brightness  error  and  a  penalty  term  such  as  a  measure 
of  departure  from  smoothness.  These  yield  second-order  partial  differential 
equations  whose  discrete  approximation  on  a  regular  grid  can  be  conveniently 
solved  by  classic  iterative  techniques  from  numerical  analysis.  Severed  of  these 
methods,  however,  compute  surface  orientation,  not  height,  and  do  not  ensure 
that  the  resulting  gradient  field  is  integrable  [Ikeuchi  &  Horn  81]  [Brooks  & 
Horn  85].  One  thus  has,  as  a  second  step,  to  find  a  surface  whose  gradient 
comes  closest  to  the  estimated  gradient  field  in  a  least-squares  sense  (see 
[Ikeuchi  84],  chapter  11  in  [Horn  86],  and  [Horn  &  Brooks  86]). 

The  two  steps  can  be  combined,  and  the  accuracy  of  the  estimated  surface 
shape  improved  considerably,  by  alternately  taking  one  step  of  the  iteration  for 
recovering  surface  orientation  from  brightness,  and  one  step  of  the  iteration 
that  recovers  the  surface  that  best  fits  the  current  estimate  of  the  surface 
gradient.  This  idea  can  be  formalized  by  setting  up  a  variational  problem 
involving  both  the  surface  height  above  a  reference  plane  and  the  first  partial 
derivatives  thereof.  The  resulting  set  of  three  coupled  Euler  equations  can  be 
discretized  and  solved  much  as  the  two  coupled  equations  are  in  the  simpler 
methods  that  only  recover  surface  orientation. 

Such  an  iterative  scheme  for  recovering  shape  from  shading  has  been 
developed  and  implemented.  The  new  scheme  recovers  height  and  gradient  at 
the  same  time.  Linearization  of  the  reflectance  map  about  the  local  average 
surface  orientation  greatly  improves  the  performance  of  the  new  algorithm 
and  could  be  used  to  improve  the  performance  of  existing  iterative  shape- 
from-shading  algorithms.  The  new  algorithm  has  been  successfully  applied  to 
complex  wrinkled  surfaces;  even  surfaces  with  discontinuities  in  the  gradient. 


9.  Acknowledgements 

Michael  Caplinger  kindly  provided  pointers  to  relevant  publications  on  pho- 
toclinometry,  in  particular  the  Ph.D.  thesis  of  R.L.  Kirk  at  the  California 
Institute  of  Technology.  Michael  Brooks  drew  my  attention  to  the  fact  that 
existing  iterative  schemes  “walk  away”  from  the  correct  solution  and  that 
regularization  does  not  “pick  one  of  the  infinite  number  of  solutions”  of  an 
ill-posed  problem,  as  sometimes  mistakenly  stated.  Bradford  Washburn  sup¬ 
plied  a  detailed  topographic  map  of  the  Mt.  Washington  area  of  the  White 
Mountains  in  New  Hampshire  from  which  digital  terrain  models  used  in  the 
experiments  here  were  interpolated.  Joel  Moses  found  a  way  for  the  Electrical 
Engineering  and  Computer  Science  Department  to  help  cover  the  cost  of  the 
author’s  travel  to  the  photoclinometry  workshop  held  at  Arizona  State  Uni¬ 
versity,  February  12-13,  1989.  This  workshop  provided  much  of  the  impetus 


56 


Height  and  Gradient  from  Shading 


for  the  work  reported  here.  Michael  Brooks,  Eric  Grimson  and  Teresa  Ehling 

provided  additional  helpful  comments. 

10.  References 

Abdou,  I.E.  $£  K.Y.  Wong  (1982)  “Analysis  of  Linear  Interpolation  Schemes  for  Bi- 
Level  Image  Applications,”  IBM  Journal  of  Research  and  Development,  Vol.  26, 
No.  6,  pp.  667-686,  November  (see  Appendix). 

Bernstein,  R.  (1976)  “Digital  Image  Processing  of  Earth  Observation  Sensor  Data,” 
IBM  Journal  of  Research  and  Development ,  pp.  40-57,  January  (see  Appendix). 

Blake,  A.,  A.  Zisserman  &  G.  Knowles  (1985)  “Surface  Descriptions  from  Stereo 
and  Shading,”  Image  &  Vision  Computing ,  Vol.  3,  No.  4,  pp.  183-191.  Also  in 
(1989)  Shape  from  Shading ,  Horn,  B.K.P.  &  M.J.  Brooks  (eds.).  MIT  Press, 
Cambridge,  MA. 

Bonner,  W.J.  &  R.A.  Schmall  (1973)  “A  Photometric  Technique  for  Determining 
Planetary  Slopes  from  Orbital  Photographs,”  U.S.  Geological  Survey  Profes¬ 
sional  Paper  812-A,  pp.  1-16. 

Brandt,  A.  (1977)  “Multi-level  Adaptive  Solutions  to  Boundary-Value  Problems,” 
Mathematics  of  Computation,  Vol.  31,  No.  138,  April,  pp.  333-390. 

Brandt,  A.  (1980)  “Stages  in  Developing  Multigrid  Solutions,”  in  Numerical  Meth¬ 
ods  for  Engineering,  Absi,  E.,  R.  Glowinski,  P.  Lascaux,  H.  Veysseyre  (eds.), 
Dunod,  Paris,  pp.  23-44. 

Brandt,  A.  (1984)  Multigrid  Techniques:  1984  Guide  with  Applications  to  Fluid  Dy¬ 
namics,  Monograph  available  as  GMD-Studie  No.  85,  from  GMD-F1T,  Post- 
fach  1240,  D-2505,  St.  Augustin  1,  West  Germany. 

Brandt,  A.  &  N.  Dinar  (1979)  “Multigrid  solutions  of  elliptic  flow  problems,”  in 
Parter,  S.V.  (ed.)  Numerical  Methods  for  PDE,  Academic  Press,  New  York, 
NY. 

Brooks,  M.J.  (1983)  “Two  Results  Concerning  Amoiguity  in  Shape  from  Shading,” 
Proceedings  of  the  National  Conference  on  Artificial  Intelligence,  Washington, 
D.C.,  August  22-26,  pp.  36-39. 

Brooks,  M.J.  (1985)  Personal  communication. 

Brooks,  M.J.  &  B.K.P.  Horn  (1985)  “Shape  and  Source  from  Shading,”  Proceedings 
of  the  International  Joint  Conference  on  Artificial  Intelligence ,  Los  Angeles, 
CA,  August  18-23,  pp.  932-936.  Also  in  (1989)  Shape  from  Shading,  Horn, 
B.K.P.  h  M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 

Bruss,  A.R.  (1982)  “The  Eikonal  Equation:  Some  Results  Applicable  to  Computer 
Vision,”  Journal  of  Mathematical  Physics,  Vol.  23,  No.  5,  pp.  890-896,  May. 
Also  in  (1989)  Shape  from  Shading,  Horn,  B.K.P.  &  M.J.  Brooks  (eds.),  MIT 
Press,  Cambridge,  MA. 


10.  References 


57 


Courant,  R.  &  D.  Hilbert  (1962)  Methods  of  Mathematical  Physics,  Volume  I,  Wiley, 
New  York,  NY. 

Davis,  P.A.  &  A.S.  McEwen  (1984)  “Photoclinometry:  Analysis  of  Inherent  Errors 
and  Implications  for  Topographic  Measurement,”  Lunar  and  Planetary  Science 
Conference  XV,  12-16  March,  pp.  194-195. 

Davis,  P.A.  &  L.A.  Soderblom  (1983)  “Rapid  Extraction  of  Relative  Topography 
from  Viking  Orbiter  Images:  II.  Application  to  Irregular  Topographic  Fea¬ 
tures,”  in  Reports  on  Planetary  Geology  Program  (1983),  NASA  Technical 
Memorandum  86246,  pp.  29-30. 

Davis,  P.A.  &  L.A.  Soderblom  (1984)  “Modeling  Crater  Topography  and  Albedo 
from  Monoscopic  Viking  Orbiter  Images:  I.  Methodology,”  Journal  of  Geo¬ 
physical  Research ,  Vol.  89,  No.  Bll,  October,  pp.  9449-9457. 

Davis,  P.A.,  L.A.  Soderblom  &  E.M.  Eliason  (1982)  “Rapid  Estimation  of  Martian 
Topography  from  Viking  Orbiter  Image  Photometry,”  in  Reports  on  Planetary 
Geology  Program  (1982) ,  NASA  Technical  Memorandum  85127,  pp.  331-332. 

Deift,  P.  &  Sylvester,  J.  (1981)  “Some  Remarks  on  the  Shape-from-Shading  Prob¬ 
lem  in  Computer  Vision,”  Journal  of  Mathematical  Analysis  and  Applications, 
Vol.  84,  No.  1,  pp.  235-248,  November. 

Frankot,  R.T.  &  R.  Chellappa  (1988)  “A  Method  for  Enforcing  Integrability  in 
Shape  from  Shading  Algorithms,”  IEEE  Transactions  on  Pattern  Analysis  and 
Machine  Intelligence,  Vol.  10,  No.  4,  pp.  439-451,  July.  Also  in  (1989)  Shape 
from  Shading,  Horn,  B.K.P.  &  M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 

Hapke,  B.W.  (1963)  “A  Theoretical  Photometric  Function  for  the  Lunar  Surface,” 
Journal  of  Geophysical  Research,  Vol.  68,  No.  15,  pp.  4571-4586,  August. 

Hapke,  B.W.  (1965)  “An  Improved  Theoretical  Lunar  Photometric  Function,”  As¬ 
tronomical  Journal,  Vol.  71,  No.  5,  pp.  333-339,  June. 

Hapke,  B.W.  (1981)  “Bidirectional  Reflectance  Spectroscopy:  (1)  Theory,”  Journal 
of  Geophysical  Research ,  Vol.  86,  No.  B4,  April,  pp.  3039-3054. 

Hapke,  B.W.  (1984)  “Bidirectional  Reflectance  Spectroscopy:  (3)  Correction  for 
Macroscopic  Roughness,”  Icarus,  Vol.  59,  pp.  41-59. 

Hapke,  B.W.  &  E.  Wells  (1981)  “Bidirectional  Reflectance  Spectroscopy:  (2)  Exper¬ 
iments  and  Observations,”  Journal  of  Geophysical  Research,  Vol.  86,  No.  B4, 
April,  pp.  3055-3060. 

Harris,  J.G.  (1986)  “The  Coupled  Depth/Slope  Approach  to  Surface  Reconstruc¬ 
tion,”  S.M.  Thesis,  Department  of  Electrical  Engineering  and  Computer  Sci¬ 
ence,  MIT.  Also  Technical  Report  908,  Artificial  Intelligence  Laboratory,  MIT, 
Cambridge,  MA. 

Harris,  J.G.  (1987)  “A  New  Approach  to  Surface  Reconstruction:  The  Coupled 
Depth/Slope  Model,”  Proceedings  of  the  International  Conference  on  Computer 
Vision,  London,  England,  June  8-11,  pp.  277-283. 


58 


Height  and  Gradient  from  Shading 


Horn,  B.K.P.  (1970)  “Shape  from  Shading:  A  Method  for  Obtaining  the  Shape  of  a 
Smooth  Opaque  Object  from  One  View,”  Ph.D.  Thesis,  Department  of  Elec¬ 
trical  Engineering,  MIT.  Also  Technical  Report  TR-79,  Project  MAC,  MIT, 
Cambridge,  MA.  Also  Technical  Report  TR-232,  Artificial  Intelligence  Labo¬ 
ratory,  MIT,  Cambridge,  MA. 

Horn,  B.K.P.  (1975)  “Obtaining  Shape  from  Shading  Information,”  Chapter  4  in 
The  Psychology  of  Computer  Vision ,  P.H.  Winston  (ed.),  McGraw  Hill,  New 
York,  NY,  pp.  115-155.  Also  in  (1989)  Shape  from  Shading,  Horn,  B.K.P.  & 
M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 

Horn,  B.K.P.  (1977)  “Understanding  Image  Intensities  (sic),”  Artificial  Intelligence, 
Vol.  8,  No.  2,  pp.  201-231,  April.  Also  in  (1987)  Readings  in  Computer  Vision, 
Fischler,  M.A.  &  O.  Firschein  (eds.),  Kaufmann,  pp.  45-60. 

Horn,  B.K.P.  (1981)  “Hill  Shading  and  the  Reflectance  Map,”  Proceedings  of  the 
IEEE,  Vol.  69,  No.  1,  pp.  14-47,  January.  Also  (1982)  Geo- Processing ,  Vol.  2, 
No.  1,  pp.  65-146,  October.  Also  (1979)  “Automatic  Hill-Shading  and  the  Re¬ 
flectance  Map,”  Image  Understanding  Workshop,  Palo  Alto,  CA,  April  24-25, 
pp.  79-120. 

Horn,  B.K.P.  (1984)  “Extended  Gaussian  Images,”  Proceedings  of  the  IEEE,  Vol.  72, 
No.  12,  pp.  1671-1686,  December. 

Horn,  B.K.P.  (1986)  Robot  Vision,  MIT  Press,  Cambridge,  MA  h  McGraw-Hill, 
New  York,  NY. 

Horn,  B.K.P.  (1988)  “Some  Ideas  on  Parallel  Analog  Networks  for  Machine  Vision,” 
Memo  1071,  Artificial  Intelligence  Laboratory,  MIT,  Cambridge,  MA,  Decem¬ 
ber. 

Horn,  B.K.P.  &:  B.L.  Bachman  (1978)  “Using  Synthetic  Images  to  Register  Real 
Images  with  Surface  Models,”  Communications  of  the  ACM,  Vol.  21,  No.  11, 
pp.  914-924,  November.  Also  “Registering  Real  Images  Using  Synthetic  Im¬ 
ages,”  in  Artificial  Intelligence:  An  MIT  Perspective  (Volume  II),  Winston, 
P.H.  &  R.H.  Brown  (eds.),  MIT  Press,  Cambridge,  MA,  pp.  129-160. 

Horn,  B.K.P.  &  M.J.  Brooks  (1986)  “The  Variational  Approach  to  Shape  from 
Shading,”  Computer  Vision,  Graphics  and  Image  Processing,  Vol.  33,  No.  2, 
pp.  174-208,  February.  Also  in  (1989)  Shape  from  Shading,  Horn,  B.K.P.  & 
M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 

Horn,  B.K.P.  &  M.J.  Brooks  (eds.)  (1989)  Shape  from  Shading,  MIT  Press,  Cam¬ 
bridge,  MA,  June. 

Horn,  B.K.P.  &  B.G.  Schunck  (1981)  “Determining  Optical  Flow,”  Artificial  Intel¬ 
ligence,  Vol.  17,  No.  1-3,  August  1981,  pp.  185-203.  Also  (1980)  Memo  572, 
Artificial  Intelligence  Laboratory,  MIT,  Cambridge,  MA,  April. 

Horn,  B.K.P.  &  R.W.  Sjoberg  (1979)  “Calculating  the  Reflectance  Map,”  Applied 
Optica,  Vol.  18,  No.  11,  pp.  1770-1779,  June.  Also  in  (1989)  Shape  from  Shading, 
Horn,  B.K.P.  &  M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 


10.  References 


59 


Howard,  A.D.,  K.R.  Blasius,  k  J.A.  Cutt  (1982)  “Photoclinometric  Determination  of 
the  Topography  of  the  Martian  North  Polar  Cap,”  Icarus ,  Vol.  50,  pp.  245-258. 

Ikeuchi,  K.  (1984)  “Reconstructing  a  Depth  Map  from  Intensity  Maps,”  Interna¬ 
tional  Conference  on  Pattern  Recognition ,  Montreal,  Canada,  July  30-August 
2,  pp.  736-738.  Also  (1983)  “Constructing  a  Depth  Map  from  Images,”  Memo 
744,  Artificial  Intelligence  Laboratory,  MIT,  Cambridge,  MA,  August. 

Ikeuchi,  K.  k  B.K.P.  Horn  (1981)  “Numerical  Shape  from  Shading  and  Occluding 
Boundaries,”  Artificial  Intelligence,  Vol.  17,  No.  1-3,  pp.  141-184,  August.  Also 
in  (1989)  Shape  from  Shading,  Horn,  B.K.P.  k  M.J.  Brooks  (eds.),  MIT  Press, 
Cambridge,  MA. 

Keys,  R.G.  (1981)  “Cubic  Convolution  Interpolation  for  Digital  Image  Processing,” 
IEEE  Transactions  on  Acoustics,  Speech  and  Signal  Processing,  Vol.  29,  No.  6, 
pp.  1153-1160,  December. 

Kirk,  R.L.  (1984)  “A  Finite-Element  Approach  to  Two-Dimensional  Photoclinome- 
try,”  brief  abstract  in  Bulletin  of  the  American  Astronomical  Society,  Vol.  16, 
No.  3,  pg.  709. 

Kirk,  R.L.  (1987)  “A  Fast  Finite-Element  Algorithm  for  Two-Dimensional  Photo- 
clinometry,”  Part  III  of  Ph.D.  Thesis,  Division  of  Geological  and  Planetary 
Sciences,  California  Institute  of  Technology,  Pasadena,  CA. 

Koenderink,  J.J.  k  A.J.  van  Doom  (1980)  “Photometric  Invariants  Related  to  Solid 
Shape,”  Optica  Acta ,  Vol.  27,  No.  7,  pp.  981-996.  Also  in  (1989)  Shape  from 
Shading,  Horn,  B.K.P.  k  M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 

Lambiotte,  J.J.  k  G.R.  Taylor  (1967)  “A  Photometric  Technique  for  Deriving  Slopes 
from  Lunar  Orbiter  Photography,”  Proceedings  of  the  Conference  on  the  Use  of 
Space  Systems  for  Planetary  Geology  and  Geophysics,  Boston,  MA,  May  25-27. 

Lee,  C.-H.  k  A.  Rosenfeld  (1985)  “Improved  Methods  of  Estimating  Shape  from 
Shading  using  the  Light  Source  Coordinate  System,”  Artificial  Intelligence, 
Vol.  26,  No.  2,  pp.  125-143.  Also  in  (1989)  Shape  from  Shading,  Horn,  B.K.P. 
k  M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 

Lee,  D.  (1988)  “Algorithms  for  Shape  from  Shading  and  Occluding  Boundaries,” 
Proceedings  of  the  IEEE  Conference  on  Computer  Vision  and  Pattern  Recogni¬ 
tion,  June  5-9,  Ann  Arbor,  MI,  pp.  478-485.  Also  in  (1989)  Shape  from  Shading, 
Horn,  B.K  r  k  H..T.  Brooks  (eds  ),  MIT  press,  Cambridge,  MA. 

Lucchitta,  B.K.  k  N.A.  Gambell  (1970)  “Evaluation  of  Photoclinometric  Profile 
Determination,”  in  Analysis  of  Apollo  8  Photographs  and  Visual  Observations, 
NASA  SP-201,  National  Aeronautics  and  Space  Administration,  pp.  51-59. 

Malik,  J.  k  D.  Maydan  (1989)  “Recovering  Three  Dimensional  Shape  from  a  Single 
Image  of  Curved  Objects,”  IEEE  Transactions  on  Pattern  Analysis  and  Ma¬ 
chine  Intelligence,  Vol.  11,  No.  6,  pp.  555-566,  June.  Also  in  (1989)  Shape  from 
Shading,  Horn,  B.K.P.  k  M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 


60 


Height  and  Gradient  from  Shading 


Malin,  M.C.  k  G.E.  Danielson  (1984)  “Topography  on  Ganymede  Derived  from 
Photoclinometry,”  in  Reports  on  Planetary  Geology  Program  (1983),  NASA 
Technical  Memorandum  86246,  pp.  29-30. 

McEwen,  A.S.  (1985)  “Topography  and  Albedo  of  Ius  Chasma,  Mars,”  Lunar  and 
Planetary  Science  Conference  XVI,  11-15  March,  pp.  528-529. 

Mingolla,  E.  k  J.T.  Todd  (1986)  “Perception  of  Solid  Shape  from  Shading,”  Bio¬ 
logical  Cybernetics,  Vol.  53,  pp.  137-151.  Also  in  (1989)  Shape  from  Shading, 
Horn,  B.K.P.  k  M.J.  Brooks  (eds.),  MIT  Press,  Cambridge,  MA. 

Minnaert,  M.  (1961)  “Photometry  of  the  Moon,”  Chapter  6  in  Volume  3  of  Planets 
and  Satellites:  The  Solar  System,  Kuiper,  G.P.  k  B.M.  Middlehurst  (eds.), 
University  of  Chicago  Press,  Chicago,  IL,  pp.  213-248. 

Passey,  Q.R.  k  E.M.  Shoemaker  (1982)  “Craters  and  Basins  on  Ganymede  and 
Callisto:  Morphological  Indicators  of  Crustal  Evolution,”  in  Morrison,  D.  (ed.) 
Satellites  of  Jupiter,  University  of  Arizona  Press,  Tucson,  pp.  379-434. 

Pentland,  A.P.  (1984)  “Local  Shading  Analysis,”  IEEE  Transactions  on  Pattern 
Analysis  and  Machine  Intelligence ,  Vol.  6,  No.  2,  pp.  170-187,  March.  Also  in 
(1989)  Shape  from  Shading,  Horn,  B.K.P.  k  M.J.  Brooks  feds.),  MIT  Tress, 
Cambridge,  MA. 

Pentland,  A.P.  (1988)  “Shape  Information  from  Shading:  A  Theory  about  Human 
Perception,”  Technical  Report  103,  Vision  Sciences,  MIT  Media  Laboratory, 
MIT,  Cambridge,  MA,  May. 

Rifman,  S.S.&  D.M.  McKinnon  (1974)  “Evaluation  of  Digital  Correction  Techni¬ 
ques — for  ERTS  images,”  Report  Number  E74-10792,  TRW  Systems  Group, 
July  (see  Chapter  4).  Also  Final  Report  TRW  20634-6003-TU-00,  NASA  God¬ 
dard  Space  Flight  Center. 

Rindfleisch,  T.  (1966)  “Photometric  Method  for  Lunar  Topography,”  Photogrammet- 
ric  Engineering,  Vol.  32,  No.  2,  pp.  262-277,  March.  Also  (1965)  “A  Photomet¬ 
ric  Method  for  Deriving  Lunar  Topographic  Information,”  Technical  Report  32- 
786,  Jet  Propulsion  Laboratory,  California  Institute  of  Technology,  Pasadena, 
CA,  September. 

Rowan,  L.C.,  J.F.  McCauley  k  E.A.  Holm  (1971)  “Lunar  Terrain  Mapping  and 
Relative  Roughness  Analysis,”  U.S.  Geological  Survey  Professional  Paper  599- 
G,  pp.  1-32. 

Saxberg,  B.V.H.  (1988)  “A  Modern  Differential  Geometric  Approach  to  Shape  from 
Shading,”  Ph.D.  Thesis,  Department  of  Electrical  Engineering  and  Computer 
Science,  MIT,  Cambridge,  MA. 

Simchony,  T.,  R.  Chellappa  k  M.  Shao  (1989)  “Direct  Analytical  Methods  for  Solv¬ 
ing  Poisson  Equations  in  Computer  Vision  Problems,”  unpublished  report,  Uni¬ 
versity  of  Southern  California. 

Sjoberg,  R.W.  k  B.K.P.  Horn  (1983)  “Atmospheric  Effects  in  Satellite  Imaging  of 
Mountainous  Terrain,”  Applied  Optics ,  Vol.  22,  No.  11,  pp.  1702-1716,  June. 


10.  References 


61 


Squyres,  S.W.  (1981)  “The  Topography  of  Ganymede’s  Grooved  Terrain,”  Icarus, 
Vol.  46,  pp.  156-168. 

Strat,  T.  (1979)  “A  Numerical  Method  for  Shape  from  Shading  for  a  Single  Im¬ 
age,”  S.M.  Thesis,  Department  of  Electrical  Engineering  and  Computer  Sci¬ 
ence,  MIT,  Cambridge,  MA. 

Sugihara,  K.  (1986)  Machine  Interpretation  of  Line  Drawings ,  Chapter  10,  MIT 
Press,  Cambridge,  MA. 

Szeliski,  R.  &  B.K.P.  Horn  (1989)  “An  Impossible  Shaded  Image,”  unpublished. 

Terzopoulos,  D.  (1983)  “Multilevel  Computational  Processes  for  Visual  Surface 
Reconstruction,”  Computer  Vision,  Graphics  and  Image  Processing ,  Vol.  24, 
pp.  52-96.  Also  (1982)  “Multi-level  Reconstruction  of  Visual  Surfaces:  Vari¬ 
ational  Principles  and  Finite  Element  Representation,  Memo  671,  Artificial 
Intelligence  Laboratory,  MIT,  Cambridge,  MA,  April. 

Terzopoulos,  D.  (1984)  “Multigrid  Relaxation  Methods  and  the  Analysis  of  Light¬ 
ness,  Shading,  and  Flow,”  Memo  803,  Artificial  Intelligence  Laboratory,  MIT, 
Cambridge,  MA,  October.  Also  chapter  10  in  Image  Understanding  84,  Ull- 
man,  S.  &:  W.  Richards  (eds.),  Ablex  Publishing  Corporation,  Norwood,  NJ, 
pp.  225-262. 

Tyler,  G.L.,  R.A.  Simpson  &  H.J.  Moore  (1971)  “Lunar  Slope  Distributions:  Com¬ 
parison  of  Bi-Static  Radar  and  Photographic  Results,”  Journal  of  Geophysical 
Research,  Vol.  76,  No.  11,  pp.  2790-2795. 

Watson,  K.  (1968)  “Photoclinometry  from  Spacecraft  Images,”  U.S.  Geological  Sur¬ 
vey  Professional  Paper  599-B,  pp.  1-10. 

Wildey,  R.L.  (1975)  “Generalized  Photoclinometry  for  Mariner  9,”  Icarus,  Vol.  25, 
pp.  613-626. 

Wildey,  R.L.  (1984)  “Topography  from  Single  Radar  Images,”  Science ,  Vol.  224, 
pp.  153-156,  April. 

Wildey,  R.L.  (1986)  “Radarclinometry  for  the  Venus  Radar  Mapper,”  Photogram- 
metric  Engineering  and  Remote  Sensing ,  Vol.  52,  No.  1,  pp.  41-50,  January. 
Also  in  (1989)  Shape  from  Shading,  Horn,  B.K.P.  &:  M.J.  Brooks  (eds.),  MIT 
Press,  Cambridge,  MA. 

Wilhelms,  D.E.  (1964)  “A  Photometric  Technique  for  Measurement  of  Lunar  Slopes,” 
in  Astrogeological  Studies  Annual  Progress  Report,  Part  D:  Studies  for  Space 
Flight  Program.  U.S.  Geological  Survey  Open-File  Report,  pp.  1-12,  May,  NASA 
Catalog  Number  N66  35597. 

Wilson,  L.,  M.A.  Brown,  E.M.  Parmentier  &  J.W.  Head  (1984)  “Theoretical  Aspects 
of  Photoclinometry  Terrain  Profiling  on  the  Galilean  Satellites,”  in  Reports 
on  Planetary  Geology  Program  (1983),  NASA  Technical  Memorandum  86246, 
pp.  27-28. 


A 


62  Height  and  Gradient  from  Shading 

Wilson,  L.,  J.S.  Hampton,  Sc  H.C.  Balen  (1985)  “Photoclinometry  of  Terrestrial 
Sc  Planetary  Surfaces,”  Lunar  and  Planetary  Science  Conference  XVI,  11-15 
March,  pp.  912-913. 

Woodham,  R.J.  (1977)  “A  Cooperative  Algorithm  for  Determining  Surface  Orien¬ 
tation  from  a  Single  View,”  International  Joint  Conference  on  Artificial  Intel¬ 
ligence,  Cambridge,  MA,  August  22-25,  pp.  635-641. 

Woodham,  R.J.  (1978)  “Photometric  Stereo:  A  Reflectance  Map  Technique  for  De¬ 
termining  Surface  Orientation  from  a  Single  View,”  Image  Understanding  Sys¬ 
tems  &  Industrial  Applications,  Proceedings  of  the  Society  of  Photo-Optical 
Instrumentation  Engineers,  Vol.  155,  pp.  136-143. 

Woodham,  R.J.  (1979)  “Analyzing  Curved  Surfaces  using  Reflectance  Map  Tech¬ 
niques,”  in  Artificial  Intelligence:  An  MIT  Perspective  (Volume  II),  Winston, 
P.H.  Sc  R.H.  Brown  (eds.),  MIT  Press,  Cambridge,  MA,  pp.  161-184. 

Woodham,  R.J.  (1980)  “Photometric  Method  for  Determining  Surface  Orientation 
from  Multiple  Images,”  Optical  Engineering,  Vol.  19,  No.  1,  January-February, 
pp.  139-144.  Also  in  (1989)  Shape  from  Shading,  Horn,  B.K.P.  Sc  M.J.  Brooks 
(eds.),  MIT  Press,  Cambridge,  MA. 

Woodham,  R.J.  (1989)  “Determining  Surface  Curvature  with  Photometric  Stereo,” 
Proceedings  IEEE  Conference  on  Robotics  and  Automation,  14- 19  May. 


